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ABSTRACT 

A new formalism is derived for the analysis and exact reconstruction of band-limited 
signals on the sphere wit h directional wavelets. It rep resents an ev olution of the wavelet 
formalism developed by lAntoine fc Vanderghevnslj (1999) and IWiaux et al.1 (|2005D . 
The translations of the wavelets at any point on the sphere and their proper rotations 
are still defined through the continuous three-dimensional rotations. The dilations 
of the wavelets are directly defined in harmonic space through a new kernel dila- 
tion, which is a modification of an existing harmonic dilation. A family of factorized 
steerable functions with compact harmonic support which arc suitable for this ker- 
nel dilation is firstly identified. A scale discretized wavelet formalism is then derived, 
relying on this dilation. The discrete nature of the analysis scales allows the exact 
reconstruction of band-limited signals. A corresponding exact multi-resolution algo- 
rithm is finally described and an implementation is tested. The formalism is of interest 
notably for the denoising or the deconvolution of signals on the sphere with a sparse 
expansion in wavelets. In astrophysics, it finds a particular application for the identi- 
fication of localized directional features in the cosmic microwave background (CMB) 
data, such as the imprint of topological defects, in particular cosmic strings, and for 
their reconstruction after separation from the other signal components. 

Key words: methods: data analysis, techniques: image processing, cosmology: cosmic 
microwave background 



■ 1 INTRODUCTION 

, Very generically, the scale-space analysis of a signal with 

■ wavelets on a given manifold defines wavelet coefficients 
' which characterize the signal around each point of the 

manifold and at va rious scales ('Mallat' 'l998'; ' Antoine et al.l 
12004 : Antoine fc Vandergheynst.2007 ). Wavel et techniques 
find n umerous applications in astrophysics (|Starck et al.l 
l2006al ). It commonly concerns the analysis of data dis- 
tributed on the real line of time, or images on the plane. 
But other experiments also acquire data in all directions 
of the sky. This is notably the case of observations of 
the cosmic microwave background (CMB) radiation, such 
as the current Wilkinson Microwave Anisotropy Probe 
(WMAP) satellite experiment, or the forthcoming Planck 
Surveyor satellite experiment. Sky surveys such as the 
NRAO Very Large Array Sky Survey (NVSS) also map 
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data on a large fraction the celestial sphere. The scale- 
space analysis of such data sets requires wavelet techniques 
on the sphere. Various wavelet formalisms have been pro- 
posed to date ( Holschncidcr 1996; Frcedon & Windhcuseij 
19961 : iFreeden et al., ,1998.: Antoine fc Vanderghevnst .19981 



19991 : iNarcowich et al.1 120051: iMcEwen et al l bood) The 
formalism originated by Antoine fc Vanderghevnst (|l999l ) 
in a gro up-theoretic contex t triggered various devel 



opments (lAntoine et al 



I2OO3I : iBogdanova et al 



I2OO2I : iDemanet fc VanderghevnstI 



2005|'). and was r econs idered in a 



more practical context by Iwiaux et all l|2005l l . This ap- 
proach notably found a recent and very interesting appli- 
cation in the a nalysis of the CMB data, as reviewed by 
iMcEwen et al.l (|2007bl ). 

In particular, the denoising or the deconvolution of data 
represents a large field of application of wavelet techniques. 
Experimental data sets are indeed always affected by var- 
ious noise sources, notably related to the instrumentation. 
The data can also be blurred by experimental beams asso- 
ciated with the instrumentation. Signals detected can also 
originate from different physical sources, which need to be 
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separated. In that component separation perspective, each 
component can in turn be technically understood as a signal, 
while the other components are seen as noise. As an exam- 
ple, observed CMB data represent a superposition of the 
CMB signal itself with instrumental noise and foreground 
emissions, blurred by the experimental beam at each detec- 
tion frequency. The denoising and the deconvolution of the 
signal, and the separation of its astrophysical components is 
of major interest for astrophysics and cosmology. 

Signals with features defined at specific positions and 
scales typically have a sparse expansion in terms of wavelets. 
For such signals, denoising and deconvolution algorithms 
are generically much more efficient when applied to the 
wavelet coefficients (Mallat 1998; Daubcchics 2004). How- 
ever this requires a scheme allowing the exact reconstruc- 
tion of the signals analyzed from their wavelet coefficients. 
Moreover, localized characteristics can be elongated, in 
which case directional wavelets are essential. The iden- 
tification and the reconstruction of localized directional 
features in CMB data represents a very interesting ap- 
plication of such a framework. It typically concerns the 
imprint of topological defects such as textures or cos- 
mic strings ^Kaiser & Stebbins"l984VTurok & Sporgel"l990|; 
IVilenkin fc Shellard.1994. : .Hind marsh fc Kibble 1995). This 
application can be recast in a component separation ap- 
proach where all continuous and typically Gaussian emis- 
sions are seen as noise, in contrast with localized directional 
features. Let us also emphasize that such a framework can 
have many applications well beyond astrophysics, from geo- 
physics to biomedical imaging, or computer vision. 

At present, the simultaneous combination of the prop- 
erties of exact reconstruction and directionality is lacking 
in the existing wavelet formalisms on the sphere. It has 
only been considered for the wavelet analysis of signals on 
the p lane (|Simoncelli et al.lll992l : IVanderghevnst fc GobbersI 
|2002| ). The primary aim of the present work resides in the 
development of a new scale discretized wavelet formalism for 
the analysis and exact reconstruction of band-limited signals 
on the sphere with directional wavelets. As a by-product, a 
new continuous wavelet formalism is also obtained, which 
allows the analysis of signals with a new family of wavelets 
relative to existing formalisms. But the continuous range of 
scales required for the analysis prevents exact reconstruction 
in practice, for which the scale discretized wavelet formalism 
proposed is essential. 

The remainder of this paper is organized as follows. In 
Section[2l we present an existing scheme for the definition of 
a continuous wavelet formalism on the sphere from a generic 
dilation operation. We consider directional and axisymmet- 
ric wavelets and discuss the cases of the stereographic and 
harmonic dilations. In Section [31 we propose a new kernel 
dilation. A corresponding family of factorized steerable func- 
tions with compact harmonic support is identified. We show 
that localization and directionality properties of such func- 
tions can be controlled through kernel dilation. In Section 
|4l we derive a new continuous wavelet formalism from the 
kernel dilation with continuous scales. We then derive a new 
scale discretized wavelet formalism that allows the exact re- 
construction of band-limited signals in practice. We design 
explicitly an example wavelet. We finally recast the scale 
discretized wavelet formalism in an invertible filter bank ap- 
proach. In Section O we describe an exact algorithm ac- 



counting for the multi-resolution properties of the formal- 
ism. The memory and computation time requirements are 
discussed and an implementation is tested. In Section [S] we 
discuss the application of the formalism to the detection of 
cosmic strings through the denoising of full-sky CMB data. 
We finally conclude in Section [T] 



2 WAVELETS FROM A GENERIC DILATION 

In this section we discuss an existing scheme for the defini- 
tion of a continuous wavelet formalism on the sphere from a 
generic dilation operation. We consider directional and ax- 
isymmetric wavelets explicitly. We finally discuss in detail 
the stereographic and harmonic dilations. 



2.1 Directional wavelets 

In the continuous fram ewor k develope d by 

lAntoine fc Vanderghev"ns^ (|l999l ') and IWiaux et al] (|2005h . 

the wavelet analysis of a signal on the sphere, i.e. the 
unit sphere S^, defines wavelet coefficients through the 
correlation of the signal with dilated versions of a local 
analysis function. Theoretically, the signal can be recovered 
explicitly from its wavelet coefficients provided that the lo- 
cal analysis functions satisfies some admissibility condition, 
raising it to the rank of a wavelet. 

The real and harmonic structures of S^ are summa- 
rized concisely as follows. We consider a three-dimensional 
Cartesian coordinate system {o,ox,oy,oz) centered on the 
sphere, and where the direction oz identifies the North pole. 
Any point co on the sphere is identified by its corresponding 
spherical coordinates {9,ip), where 6 £ [0, tt] stands for the 
co-latitude, or polar angle, and ip G [0, 2-n) for the longitude, 
or azimuthal angle. Let the continuous signal Fiua) and the 
local analysis function $(a;) be square-integrable functions 
on the sphere: F, ^ £ L'^(S'^, df2), with the invariant mea- 
sure AO, — dcosSd^s. The spherical harmonics form an or- 
thonormal basis for the decomposition of square-integrable 
functions. They are explicitly given in a factorized form in 
terms of the associated Legendre polynomials P[^{coa6) and 
the complex exponentials e''"^^ as 



Ylr. 



2l + l{l~ my. 



-,1/2 



47r {l + my. 



Pr (cos e) , 



(1) 



with Z g N, m £ Z, and \m\ < I (|Abramowitz fc StegunI 
1 19651 : IVarshalovich et aLlll989^ . The index I represents an 
overall frequency on the sphere. The absolute value \m\ rep- 
resents the frequency associated with the azimuthal vari- 
able if. Any function G £ L^(S^,dr2) is thus uniquely 
given as a linear combination of scalar spherical harmonics: 
G{uj) = Y.ien'l^imi^iGimYimiuj). This combination defines 
the inverse spherical harmonic transform on S^. The cor- 
responding spherical harmonic coefficients are given by the 
projection Gim = {Yim\G), with \m\ ^ I, where the bracket 
{F2\Fi) — /g2 dQF2{uj)Fi{uj) generically denotes the scalar 
product for F\,F2 G L^(S^,df2). This projection defines the 
direct spherical harmonic transform on S^. 

Continuous affine transformations such as translations, 
rotations, and dilations are applied to the analysis function. 
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The continuous translations by luq — {do,ipo) G and ro- 
tations by X G [0, Stt) are defined by the three Euler an- 
gles defining an element p = (1^0,^0, x) of the group of ro- 
tations in three dimensions SO(3). The operator R{ujo) in 
L'^(S'^,dfl) for the translation of amplitude loq = {9Q,ipo) of 
a function G reads as 

(to) = [R (luo) G] (c^) = G , (2) 

where R^g{9,Lp) — [R^^Rg^]{9,if) is defined by the three- 
dimensional rotation matrices Rg^ and R^^ , acting on the 
Cartesian coordinates {x,y,z) associated with ui — {0,ip). 
The rotation operator -R^(x) in L (S ,df2) for the rotation 
of the function G around itself, by an angle x G [0, 2n), is 
given as 

(cu) = [r' (x) g] (u;) = G (r^'u;) , (3) 

where R^{0, tp) = (6,(p + x) ^Iso follows from the action of 
the three-dimensional rotation matrix on the Cartesian 
coordinates {x,y,z) associated with ui — {0,Lp). The opera- 
tor incorporating both the translations and rotations simply 
reads as Rip) = R{u)q)R\x) and Gp{u) = [R{p)G\{u) = 
G{Rp^Lo), with Rp — RuiqRx- The continuous dilations af- 
fect by definition the continuous scale of the function. The 
notion of scale may a priori be defined both in real or in har- 
monic space on S'^. In the remainder of the present subsec- 
tion we simply denote the dilated function as Ga{u}), where 
a € R+ stands for a continuous dilation factor. We explicitly 
discuss two possible definitions of dilations in Subsections 
[23] and EH 

The analysis of the signal F with an analysis function 
5' defines wavelet coefficients through the directional cor- 
relation of F with the dilated functions "Ja, i.e. the scalar 
products 

Wi{p,a) = {-i'p,a\F). (4) 

At each scale a, the wavelet coefficients W^{p,a) there- 
fore identify a square-integrable function on the rotation 
group in three dimensions SO (3). They characterize the sig- 
nal around each point too, and in each orientation x- This 
defines the scale-space nature of the wavelet decomposition 
on the sphere. 

The real and harmonic structures of the rotation group 
in three dimensions 80(3) are summarized concisely as fol- 
lows. As discussed, any rotation p on SO(3) is given in terms 
of the three Euler angles p — {ip,9,x), with 9 G [0,7r], and 
X G [0, 27r). Let H{p) be a square-integrable function 
on SO(3): H G L^(SO(3), dp), with the invariant measure 
dp = d(/pdcosMx. The Wigner D- functions are the matrix 
elements of the irreducible unitary representations of weight 
I of the group in L^(S0(3), dp). By the Peter- Weyl theorem 
on compact groups, the matrix elements Dl^„ also form an 
orthogonal basis in L^(SO(3), dp). They are explicitly given 
in a factorized form in terms of the real Wigner d-functions 
dlnniS) and the complex exponentials, e""""^ and e~^"^, as 

Dl„i^,9,x)^e~'"^''dl„{9)e-^"^, (5) 

with ^ G N, TO, n G Z, and |m |, \n\ ^ I (|Varshalovich et al.l 
Il989l : iBrink fc Satchleij|l993l ). Again, I represents an over- 
all frequency on SO(3), and \m\ and |n| the frequen- 
cies associated with the variables ip and Xi respectively. 
Any function H G L^(S0(3),dp) is thus uniquely given 



as a linear combination of Wigner D-functions: H{p) = 
EieN(2^ + l)/8vr2E|,„|,|„Ki-ffmn-Dmn(p). This combination 
defines the inverse Wigner D- function transform on SO (3). 
The corresponding Wigner D-function coefficients are given 
by the projection i/^„ = J^^^^^ dp Dl„^{p)H (p) . This pro- 
jection defines the direct Wigner D-function transform on 
SO(3). 

At each scale, the direct Wigner D-function transform 
of the wavelet coefficients is given as the pointwise product 
of the spherical harmonic coefficients of the signal and the 
wavelet: 

mtja)^^{^)lF^. (6) 

Indeed, the orthonormality of scalar spherical har- 
monics implies the Plancherel relation (-f2|^'i) = 
E,6nE|™|^, (K)L(K),„ for Fi,F2 G L'HS^,dn), and the 
action of the operator R{p) on G G L^(S^,dr2) reads in 
terms of its spherical harmonic coefficients as {Gp)i^ = 

E|„|^; Dlnn{p)Gln- 

The reconstruction of a signal F from its wavelet coef- 
ficients with an analysis function 'I' is given as 

F{lo)= f dti{a) f dpWi {p,a)[R{p)L^'fa]{Lo). 

JR*_^ iS0(3) 

(7) 

In this relation, the scale integration measure dp{a) is part 
of the definition of the dilation operation itself (see Sub- 
sections [2]3] and [2]4|. The operator in L^(S^,dn) is de- 
fined by its action on the spherical harmonic coefficients of a 
function G: LmGim = Gim./C\i. The reconstruction formula 
holds if and only if the analysis function satisfy the following 
admissibility condition for alH G N: 

° < ^* = 2rTT ^ /- '^^^"'^ < ~- (8) 

In this case the analysis function 'I' is by definition raised to 
the rank of a wavelet. From relation ((Sjl, the admissibility 
condition intuitively requires that the whole wavelet family 
{5'a(ti')}, for a G R+, covers each frequency index I with a 
finite and non-zero amplitude, hence preserving the signal 
information at each frequency. Notice that a direct connec- 
tion exists between the generic relations ([T]) and (O for the 
signal reconstruction, an d the theory of frames on the sphere 
jBogdanova et al.llioosl '). 

We generally consider band-limited signals. Any func- 
tion G G L^(S^,dn) is said to be band-limited with band 
limit B, for any B G N°, if Gim = for all I, m with I ^ B. 
Any function H G L^(SO(3),dp) is said to be band-limited 
with band hmit B, for any B G N°, if = for aU 

l,m,n with I ^ B. From relation ([6l, if the signal F or the 
wavelet ^ are band-limited on S^, then the wavelet coeffi- 
cients are automatically band-limited on S0(3), with 
the same band limit B. 

Let us already notice that the reconstruction is en- 
sured theoretically from relation (O, through the integra- 
tion on the continuous parameter p G S0(3) for trans- 
lations and rotations of the wavelet, and on the contin- 
uous dilation factor a G R+. But in practice, the re- 
construction would require the definition of exact quadra- 
ture rules for the numerical integrations. Exact quadra- 
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ture rules for integra tion of band-limited s ignals on 
exist on eg ui-angular (|Driscoll fc Healv 1 1994 ) and Gauss- 
Legendre (|Doroshkevicli et alJ l20p5a t ) pixelization s of 
{9o,<po)- HEALPix pixelizations (|G6rski et al.1 l2005l fl of 
{9o,(po) on provide approximate quadrature rules which 
can also be made very precise thanks to an iteration process. 
Pixelizations may for instance be defined on SO{S) by com- 
bining pixelizations on with an equi-angular sampling of 
X- Corresponding quadrature rules can be made exact on 
the pixelizations based on equi-angular and Gauss-Legendre 
pixelizations on S^, while those based on HEALPix pixeliza- 
tions are approximate. This extension basically relies on the 
separa ti on of the in tegration variables (|Maslen fc Rockmord 
ll997al lbl: lKostelec fc Rockmorc 2003) from relation (O. 

However, exact quadrature rules do not exist for the 
integration over scales a £ R+. In practice, this prevents 
an exact reconstruction of the signal analyzed. A scheme 
allowing an exact reconstruction requires a discretization of 
the dilation factor. A scale discretized wavelet formalism is 
proposed in Section|4l thanks to a specific choice of dilation, 
and through an integration of the dilation factor by slices in 
relation ((Tj. 



2.2 Axisymmetric wavelets 

Any general function G £ h^{S'^,dQ,) explicitly depen- 
dent on the azimuthal angle ip, is said to be directional: 
G = G{6,(p). By opposition, any function A £ L^(S^,df2) 
independent of the azimuthal angle if is said to be zonal, or 
axisymmetric: A = A{9). It only exhibits non-zero spherical 
harmonic coefficients for m — 0: Aim = AioSmO- 

In this particular case, the directional correlation of a 
signal F with A reduces to a stand ard correlation obv iously 
independent of the rotation angle x (| Wiaux et al.l200^ '). The 
analysis of F with an axisymmetric analysis function A de- 
fines wavelet coefficients through the standard correlation of 
F with the dilated functions Aa, i.e. the scalar products 



I € N: 



WK{uo,a)^{A^„a\F). 



(9) 



At each scale a, the wavelet coefficients identify a square- 
integrable function on rather than on SO (3). The spher- 
ical harmonic transform of the wavelet coefficients is still 
given as the pointwise product of the spherical harmonic 
coefficients of the signal and the wavelet: 



{Aa)ioFlr, 



(10) 



This relation simply follows from relation ((6)1 and the equal- 
ity 7?Lo(^,0) = [An/{21 + l)]i/2y;^(^). 

The reconstruction of F from its wavelet coefficients 
reads as: 

F (w) = / d^ (a) / dtjo Wa (t^o, a) [R (cjq) La A„] (uj) , 

(11) 

for any scale integration measure dfi{a), and with the op- 
erator La in L^(S^dn) defined by: lTGiq = Gio/Ca. The 
reconstruction formula holds if and only if the analysis func- 
tion satisfies the following admissibility condition for all 



http: / /healpix. jpl.nasa.gov/ 



0<C'a 



47r 
21 + 1 



dfi{a) \{Aa)i, 



< oo. 



(12) 



2.3 Stereographic dilation 

In the original set up proposed bv lAntoine fc VanderghevnstI 
( 1999,) the stereographic dilation of functions is considered, 
which is explicitly defined in real space on S^. The stere- 
ographic dilation operator D{a) on G € L'^(S^,df2), for a 
continuous dilation factor a £ R!j. , is defined in terms of the 
inverse of the corresponding stereographic dilation Da on 
points in S'^. It reads as 



Ga (uj) = [D (a) G] (ij) 

= \^^^ {a,e)G{D-^Lj) . 



(13) 



with A^/^ (a, 6*) = a-^[l+tan^ {e/2)]/[l + a-^ tan^{e/2)]. The 
dilated point is given by Da{0,ip) = (6a(0),Lp) with the lin- 
ear relation taii{9a{0) /2) — atan(6/2). The dilation oper- 
ator therefore maps the sphere without its South pole on 
itself: 9a{9) ■ 9 £ [0,-k) ^ 9a € [0,7r). This dilation operator 
is uniquely defined by the requirement of the following natu- 
ral properties. The dilation of points on S^ must be a radial 
(i.e. only affecting the radial variable 9 independently of Lp, 
and leaving invariant) and conformal (i.e. preserving the 
measure of angles in the tangent plane at each point) diffeo- 
morphism (i.e. a continuously differentiable bijection). The 
normalization by X^^^{a, 9) in (|13p is uniquely determined by 
the requirement that the dilation of functions in L^(S^, df2) 
be a unitary operator (i.e. preserving the scalar product in 
L^{S^,dQ), and specifically the norm of functions). Notice 
that the stereographic dilation operation is supported by a 
group structure for the composition law of the correspond- 
ing operator D{a). A group homomorphism also holds with 
the operation of multiplication by a on . 

In this setting, the effect of the dilation on the spheri- 
cal harmonic coefficients of the dilated function is not eas- 
ily tractable analytically. Consequently, the admissibility 
condition ((SJ is difficult to check in practice. On the con- 
trary, wavelets on the plane are well-known, and may be 
easily constructed, as the corresponding admissibility con- 
dition reduces to a zero mean condition for a function both 
integrable and square-integrable. In that context, a corre- 
spondence principle was proved (|Wiaux et al.ll2005l ). stat- 
ing that the inverse stereographic projection of a wavelet 
on the plane leads to a wavelet on the sphere. This cor- 
respondence principle notably requires the definition of a 
scale integration measure identical to the measure used 
on the plane: d^(a) = a~'^da. Notice that this measure 
naturally appears in the origin al group-theoretic context 
(|Antoine fc Vanderghevnstll 19991 ). 



2.4 Harmonic dilation 

Another possible definition of the dilation of functions may 
be considered, which is explicitly defined in harmonic space 
on S^. It was proposed in previous developments relative 
to the definition of a wavelet formalism on the sphere 
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jHolschneideill 19961 : [McE wen et al.ll2006l ). The harmonic di- 
lation is defined directly on G £ L^(S'^,dn) through a se- 
quence of prescriptions rather than in terms of the applica- 
tion of an simple operator. Firstly, an arbitrary prescription 
must be chosen to define a set of generating functions Gm{k) 
of a continuous variable k £ R+, for each m £ Z. These func- 
tions are identified to the spherical harmonic coefficients of 
G through: Gm{l) = Gim for I £ N, and \m\ ^ Secondly, 
the variable k is dilated linearly, k — I ^ k = al, just as 
would be the norm of the Fourier frequency on the plane. 
For a continuous dilation factor a £ R^^-, the spherical har- 
monic coefficients of the dilated function Ga are defined by: 

(gT),™ = (al) . (14) 

In the corresponding continuous wavelet formalism, the 
analysis function ^' must satisfy the following form of the 
admissibility condition (|8]). On the one hand \['oo = ^'o (0) = 
0, which corresponds to the requirement that ^' has a zero 
mean on the sphere: 

-!- / df7*(tj) = 0. (15) 

This zero mean is of course preserved through harmonic dila- 
tion. As the zero frequency is not supported by the wavelets, 
only signals with zero mean can be analyzed in this formal- 
ism (see relation Let us remark that wavelets on the 
sphere dilated through the stereographic dilation do not nec- 
essarily have a zero mean. On the other hand, the scale inte- 
gration measure can arbitrarily be chosen as d/i(a) = a~^da. 
This leads to a simple expression of the remaining con- 
straints for Z £ N" as 

0<G;^^E/ ^l*^(^')r<oo. (16) 

The left-hand side inequality implies < 
dk' /k' \'^rno{k')\'^ for at least one of the first two 
generating functions: mo £ {0, 1}. In other words, ei- 
ther v&o or ^'i must be non-zero on a set of non-zero 
measure on R+. The right-hand side inequality implies 
J^^dk' /k' \'^rn{k')\^ < oo for all generating functions: 
m £ Z. Hence, the generating functions must satisfy 
^m(O) = (this condition encompasses the zero mean 
condition 1)15^ in the form ^o(O) = 0) and tend to zero when 
fc' ^ oo. With this choice of scale integration measure, 
the constraints summarize to the requirement that each 
generating function satisfies a condition very similar to the 
wavelet admiss ibility condition for a n axisymmetric wavelet 
on the plane (|Antoine et al.l l2004| j^ defined by a Fourier 
transform identical to (k). Consequently, the wavelet 
admissibility condition (|16|) can be checked in practice 
and wavelets associated with the harmonic dilation can be 
designed easily. 

For continuous axisymmetric wavelets, a unique gen- 
erating function Ao{k) of a continuous variable k £ R+ is 
required. The admissibility condition p2p reduces to the fol- 
lowing expression. The analysis function A must have a zero 

The exact wavelet admissibility condition on the plane reduces 
to a zero mean condition for functions that are both integrable 
and square-integrable. 



mean and only allows the analysis of signals with zero mean. 
A unique additional condition holds independently of I: 



< Ca = 




This condition actually encompasses the zero mean condi- 
tion in the form Ao{0) — 0, and also requires that the gen- 
erating function must tend to zero when fc' — > oo. The co- 
efficients entering the reconstruction formula ([TTJ read as 
G\ = AnGA/{2l -)- 1), for « £ N°. 

2.5 Discussion 

On the one hand, the harmonic dilation lacks some of the 
important properties which hold under stereographic dila- 
tion. As the harmonic dilation does not act on points, the 
question of the corresponding properties of a radial and con- 
formal diffeomorphism make no sense. The harmonic dila- 
tion of functions is not either a unitary procedure. It does 
not preserve the scalar product in L'^(S'^, df2), or specifically 
the norm of functions. This is due to the requirement of defi- 
nition of generating functions for any function to be dilated. 
A group structure for the composition of harmonic dilations 
holds only if successive dilations of a function G are defined 
through linear dilation of the variable fc of a unique gener- 
ating function Gm(fc) for each m £ Z. The same condition 
applies for the existence of a corresponding homomorphism 
structure with the operation of multiplication by a on R!j_. 

Moreover, the harmonic dilation is explicitly defined 
in harmonic space. The evolution in real space of localiza- 
tion and directionality properties of functions on the sphere 
through harmonic dilation is therefore not known analyti- 
cally. However, in the Euclidean limit where a function is 
localized on a small portion of the sphere, this portion is as- 
similated to the tangent plane, and the stereographic and 
harmonic dilations both identify with the standard dila- 
tion in the plane. For each m £ Z, the overall frequency 
index Z £ N ^ oo identifies with the continuous variable 
k £ R+ oo, correspo nding to the norm of the Fourier 
frequency on the plane (|Holschneideij [l996t ) . So in partic- 
ular, the evolution of localization properties of functions 
through harmonic dilation is at least controlled in the Eu- 
clidean limit. 

On the other hand, the very simple action of the har- 
monic dilation in harmonic space also exhibits several ad- 
vantages relative to the stereographic dilation. Notably, the 
harmonic dilation ensures that the band limit of a wavelet 
and of the corresponding wavelet coefficients, is reduced by 
a factor a. Such a multi-resolution property is essential in 
reducing the memory and computation time requirements 
for the wavelet analysis of signals. It does not hold under 
stereographic dilation. Moreover, as already emphasized, a 
scheme allowing an exact reconstruction of signals from their 
wavelet coefficients requires a discretization of the dilation 
factor. The definition of a scale discretized wavelet formalism 
through an integration of the dilation factor a by slices in 
the continuous wavelet formalism turns out to be very nat- 
ural with a dilation defined in harmonic space, but not with 
the stereographic dilation. Indeed, one would like the dila- 
tion operation acting on scale discretized functions after the 
integration of the dilation factor by slices to be the same as 
the original dilation operation. It will become obvious that 
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this property holds for a dilation defined in harmonic space, 
but not for the stereographic dilation. 

In conclusion, no obvious definition of dilation is im- 
posed for the development of a wavelet formalism on the 
sphere. But considering our aim for a scale discretized 
wavelet formalism, as well as the essential criterion of defin- 
ing a formalism with multi-resolution properties, we will fo- 
cus on a scale discretized wavelet formalism from a dilation 
defined in harmonic space. However, for any dilation de- 
fined in harmonic space, the evolution of the localization and 
directionality properties of functions in real space through 
dilation needs to be understood and controlled. In that re- 
gard, we amend the harmonic dilation (|14|l and define a 
kernel dilation to be applied on functions which are said 
to be factorized steerable functions with compact harmonic 
support. Moreover, the kernel dilation will also render the 
transition between the continuous and scale discretized for- 
malism much simpler and more transparent than what the 
harmonic dilation can provide. 



3 KERNEL DILATION 

In this section we define the kernel dilation on factorized 
functions in harmonic space on the sphere. We consider in 
particular factorized steerable functions with compact har- 
monic support. We also study the localization and direction- 
ality properties in real space for such functions, as well as 
the controlled evolution of these properties through kernel 
dilation. 



3.1 Factorized functions and kernel dilation 

A function G e L2(S^d^7) can be defined to be a factorized 
function in harmonic space if it can be written in the form: 



Glm = Kg (0 Sir, 



(18) 



for I £ N, and |m| ^ I. The positive real kernel Kaik) G R+ 
is a generating function of a continuous variable k G R+, 
initially evaluated on integer values k = I. The directionality 
coefficients Sf^ , for / € N, and \m\ ^ I, define the directional 
split of the function. In particular, for a real function G, they 
bear the same symmetry relation as the spherical harmonic 
coefficients Gim themselves: Sj^ = {—l)™'S'i 
loss of generality one can impose 



fl_^y Without 



E 



Sj^l^ — 1, 



(19) 



for the values of I for which Sf^ is non-zero for at least one 
value of m. Hence localization properties of a function G, 
such as a measure of dispersion of angular distances around 
its central position as weighted by the function values, are 
governed by the kernel and to a lesser extent by the direc- 
tional split. Indeed, the power contained in the function G 
at each allowed value of I is fixed by the kernel only. The 
norm of G G L2(S^dn) reads as \\G\f = ^-^(0, where 

the sum runs over the values of / for which is non-zero 
for at least one value of m. However, the directional split is 
essential in defining the directionality properties measuring 
the behaviour of the function with the azimuthal variable 
</3, because of it bears the entire dependence of the spherical 
harmonic coefficients of the function in the index m. 



The kernel dilation applied to a factorized function (|18|l 
is simply defined by application of the harmonic dilation (|14|) 
to the kernel only. The directionality of the dilated function 
is defined through the same directional split as the original 
function. For a continuous dilation factor a £ R+, the dilated 
function therefore reads as: 



(Ga),„ = Kg (al) S?^ 



(20) 



Let us emphasize that the directionality coefficients 
are not affected by dilations, on the contrary of what the 
complete action of the harmonic dilation (|14|l would imply. 
The kernel and harmonic dilations strictly identify with one 
another when applied to factorized axisymmetric functions 
A, for which the directional split takes the trivial values 
SL = 5mO for / G N. 



3.2 Compact harmonic support 

Any function G G L^(S^, df2) can be said to have a compact 
harmonic support in the interval I G ([a~^i3j ,B), for any 
_B G and any real value a > 1, if 

Gim=0 for all l,m with Z ^ ( [a^^Sj , B) , (21) 

where [ij denotes the largest integer value below x G R. 
Notice that the compactness of the harmonic support of G 
can be defined as the ratio of the band limit to the width of 
its support interval. 

For a factorized function G of the form (|18|) . the com- 
pact harmonic support in the interval I G ([a~^i3j ,B) is 
ensured by the choice of a kernel with compact support in 
the interval k G (a'^B, B): 



KG{k)^0 for ki{a.-^B,B). 



(22) 



The compactness of the harmonic support of G can sim- 
ply be estimated from the compact support of the kernel as 
c(a) = a/(a — 1) G [1, oo). One has c{a) — > oo when a — > 1, 
and c(a) 1 when a oo. Typical values would be a = 2 
corresponding to a compactness c(2) = 2, or a = 1.1 leading 
to a higher compactness c(l.l) = 11. 

By a kernel dilation with a dilation factor a G R+ in 
(|20|) . the compact support of the dilated kernel Kg {ak) G 
R+ is defined in the interval k G (a~^a' -^B,a-^B). The 
compact harmonic support of the dilated function Ga it- 
self is thus defined in the corresponding interval / G 
([a^^a^^Sj , [a~^_B]), where [a;] denotes the smallest inte- 
ger value above a; G R. In particular, the compactness of the 
harmonic support of a function remains invariant through a 
kernel dilation. 



3.3 Steerable functions 

The notion of s tccrability was first introduced on the plane 
(|Freeman fc Ad clson 199ll: [Simonc clli ct al. 1992), and more 
recently defined on the sphere ( Wiaux et al.ll2005l ). By def- 
inition, G G L'^(S^,df2) is steerable if any rotation of the 
function around itself may be expressed as a linear combi- 
nation of a finite number Al of basis functions G^: 



(23) 
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The square-integrable functions kp{x) on the circle S^, with 
1 ^ m ^ M, and Al £ N", are caUed interpolation 
weights. Intuitively, steerable functions have a non-zero an- 
gular width in the azimuthal angle tp, which renders them 
sensitive to a range of directions and enables them to satisfy 
the steerability relation. This non-zero angular width natu- 
rally corresponds to an azimuthal band limit N G N'^ in the 
frequency index m associated with the azimuthal variable 
ip: 

Gim=0 for aU l,m with \m\^N. (24) 

It can actually be shown that the property of steerability 
H23p is equivalent to the existence of an azimuthal band limit 

N dm. 

On the one hand, if a function G is steerable with M 
basis functions, then the number T of values of m for which 
Gitu has a non-zero value for at least one value of I is lower 
or equal to M: M ^ T. This was firstly estab lished for 
functions on the plane ([Freeman fc AdelsonlllQQll ). and the 
proof is absolutely identical on the sphere. As a consequence, 
the function has some azimuthal band limit A'^, with T ^ 
2N - 1. 

On the other hand, if a function G has an azimuthal 
band limit A'', then it is steerable, and the number of basis 
functions can be reduced at least to M = 2N — 1. This 
second part of the equivalence can be proved by explic- 
itly deriving a steerability relation for band-limited func- 
tions with an azimuthal band limit A'^. Any band-limited 
function G can in particular be steered using M rotated 
versions G^p = R^{xp)G as basis functions, and interpola- 
tion weights given by simple translations by Xp of a unique 
square-integrable function k{x) on the circle S^: 

A/-1 

G^H= ^ fe(x-Xp)GxpH, (25) 

for specific rotation angles Xp with ^ p ^ Af — 1. One 
may choose M = 2N — 1 equally spaced rotation angles 
Xp e [0, 27r) as Xp = 2Tvp/{2N - 1), with p < 2iV - 2. 
The function fe(x) is then defined by the Fourier coefficients 
km = 1/{2N — 1) for \m\ ^ N — 1 and k,n = otherwise. 
Notice that the angles Xp s^nd the structure of the function 
fe(x) are independent of the explicit non-zero values Gi,n- 

Typically, if Gim has a non-zero value for at least one 
value of I for aU m with |m| < iV - 1, then T = 2iV - 1 
and the function is optimally steered by these M = T an- 
gles and the function fc(x) described. On the contrary, when 
values of m, with \m\ ^ TV — 1, exist for which Gim = 
for all values of I, then T < 2N — 1 and one might want to 
reduce the number A'l = 27V — 1 of basis functions. Depend- 
ing on the distribution of the T values of m for which Gi^ 
has a non-zero value for at least one value of I, the number 
of basis functions required to steer the band-limited func- 
tion may indeed be optimized to its smallest possible value 
M = T. This optimization is notably reachable for functions 
with specific distributions of the T values of m, correspond- 
ing to particular symmetries in real space. For example, a 
function G is even or odd through rotation around itself by 
X = TT if and only if Gim has non-zero values only for, re- 
spectively, even or odd values of m. This property notably 
implies that the central position of the function G identifies 
with the North pole, in the sense that its modulus \G\ is 



then always even through rotation around itself by x = 
The combination of an azimuthal band limit A*' with that 
symmetry reads as: 

Gim = for all l,m with m ^ Tiv, (26) 

with 

Tn = {-{N-1),-{N-3), {N-3),(N- 1)} . (27) 

In this particular case, T — N and one may choose M = N 
equally spaced rotation angles Xp £ [0, tt) as Xp ~ t^p/N , 
with Q ^ p ^ N — 1, and steer the function through relation 
(|25|l . The function fc(x) is defined by the Fourier coefficients 
km = 1/Af for m G Tjv and km ~ otherwise. 

In summary, the property of steerability is indeed equiv- 
alent to the existence of an azimuthal band limit in m. 
For a factorized function G of the form (|18|) . steerability 
constraints such as (|24p and (|26|) are ensured by the direc- 
tionality coefficients 5;*^, independently of the kernel. Con- 
sequently, any relation of steerability remains unchanged 
through a kernel dilation (|2Up. which by definition only af- 
fects the kernel. 

3.4 Localization control 

Let us consider the Euclidean limit where a function is lo- 
calized on a small portion of the sphere which can be assim- 
ilated to the tangent plane. As discussed, the harmonic dila- 
tion (|14p id entifies with the sta ndard dilation in the plane in 
that limit (|Holschneidedll996l ). Hence, for factorized steer- 
able functions with compact harmonic support, the kernel 
dilation (|20p certainly shares the same property if it iden- 
tifies with the harmonic dilation itself in the limit / — > oo. 
This is ensured by considering functions with directionality 
coefficients Sf^ which become independent of I in the limit 
I —> oo. Consequently, the evolution of localization proper- 
ties of functions through kernel dilation is also controlled in 
the Euclidean limit. But a much more important localization 
property holds for the kernel dilation at any frequency range 
for factorized functions with compact harmonic support. 

A typical localization property of a function G £ 
L2(S^d^) is a measure of dispersion of angular distances 
around its central position, as weighted by the function val- 
ues. The corresponding measure in harmonic space is defined 
by the dispersion of the values of / around their central posi- 
tion, as weighted by the values of the spherical harmonic co- 
efficients Gim, for each value of m. It is well-known that the 
smaller the dispersion in real space, the larger the dispersion 
in harmonic space. An optimal Dirac delta distribution on 
the sphere (5g2(tj) exhibits an infinite series in I of spherical 

harmonic coefficients: (^s^);™ = [(2^ + l)/47r]^/^5mO- On the 
contrary a spherical harmonic Yim, completely non- localized 
in real space on , by definition exhibits a unique frequency 
I. 

In particular, we need to understand the evolution of 
this localization property of a factorized steerable function 
G with compact harmonic support through the kernel di- 
lation (|20p . Let us consider for simplicity a factorized ax- 
isymmetric function A with compact harmonic support, for 
which the kernel and harmonic dilations identify with one 
another. For an initial compact harmonic support in the in- 
terval I £ ([q~^_BJ ,-B), the kernel dilation by a factor a 
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modifies the interval to I G ([a^^a^^Bj , [a^^S]). Hence 
in harmonic space the width of the harmonic support in- 
terval is multiplied by . This also measures the evolu- 
tion of the dispersion in harmonic space. In real space, one 
can intuitively consider that the corresponding dispersion of 
the values of the angular distance 6 around the North pole 
(which is the central position of any axisymmetric function) 
is multiplied by a. This intuition is actually only exact in 
the Euclidean limit I — » oo, reached when a ^ 0. But a 
weaker property holds though, on a wide class of axisym- 
metric functions on the sphere, in particular on factorized 
axisymmetric functions with compact harmonic support. It 
takes the form of the following upper bound through the 
kernel dilation by o of such an axisymmetric function A at 
a given angular distance 6 from the North pole: 

(g) K 6(^.,) " ' (28) 
1 + {0/a) 

for any inte ger fc ^ 2 and for some constant 6(71, fc) depending 
on A and k (|Narcowich et al1l2005l ). The ratio of the bounds 
at the North pole and at any fixed angular distance 6 simply 
reads as l + {6/a)^ . When a increases, this ratio gets closer to 
unity and the bound is less constraining, enabling a larger 
dispersion of the values of the angular distance 9 around 
the North pole. When a decreases, the ratio increases and 
the bound is more constraining, hence imposing a smaller 
dispersion of the values of the angular distance 6 around the 
North pole. This ensures a good behaviour in real space for 
the kernel dilation, when applied to factorized axisymmetric 
functions with compact harmonic support. 

In summary, the dispersion of angular distances around 
the central position of a function G defines a localization 
property. We have shown that the evolution of the local- 
ization of factorized axisymmetric functions with compact 
harmonic support through kernel dilation is controlled by 
the bound (|28p . For completeness, the corresponding bound 
should be analyzed for the kernel dilation of factorized steer- 
able functions with compact harmonic support, but this 
goes beyond the scope of the present work. The verifica- 
tion of more detailed localization properties in real space 
for a function designed from its spherical harmonic coeffi- 
cients requires a numerical evaluation of sampled values of 
that function. 

3.5 Directionality control 

Let us consider a typical directionality property of a function 
G € L'^(S'^,dO), such as measured by its auto-correlation 
function. The auto-correlation function of G is defined as the 
scalar product between two rotated versions of the function 
around itself by angles XiX' G [0,27r). This auto-correlation 
only depends on the difi'erence of the rotation angles Ax = 
X — x' is therefore considered in the space L^(S^,dx) 
of square-integrable functions on the circle S^: C'~'(Ax) = 
(Gx|G^'). The peakedness of the auto-correlation function 
in Ax can be considered as a measure of the directionality 
of the function: the more pea ked the auto-corre lation, the 
more directional the function ( Wiaux et ahllioosh . From the 
Plancherel relation (F^lFi) = E^^m E|™k, (?^)yK),„, 
for Fi,F2 G L^(S'^,df2), and the expression {Gx)im ~ 
e~'^"^^Gim for the action of the operator R^{x) on G, one 



gets 

C^(Ax) = ^ ^ e-""^^|G,„|^ (29) 

The value of the auto-correlation function at Ax = ob- 
viously defines the square of the norm of the function : 
G^ (0) = !|G||^ 

For a factorized function ([18]), the auto-correlation func- 
tion is strongly related to the directional split. Let us also re- 
call that in the case of a steerable function G defined through 
(|25l) . the interpolation weights depend on the values of m 
for which the spherical harmonic coefficients have non-zero 
values and on the rotation angles Xp, but not on the val- 
ues Gim themselves. This leaves enough freedom to design 
a suitable auto-correlation function and thus control the di- 
rectionality of the function. Let us also consider a compact 
harmonic support (|2ip in the interval ([q:~^_BJ ,B). We an- 
alyze the particular case where the directionality coefficients 
are independent of I for I ^ N — 1, 

SL^Sfk-i)rn for all l,m with I ^ N ~ 1, (30) 

and where A'' — 1 is lower or equal to the lowest integer value 
above the lower bound of the compact harmonic support 
interval, i.e. N — I ^ [a^^Bj -|- 1. In that limit, the auto- 
correlation reads as 

G^(Ax) = ||G||^ J2 e"™^'^!^,^!)™!'. (31) 

|m|!j JV-l 

In other words, the square of the complex norm of the direc- 
tionality coefficients identifies with the Fourier coefficients 
of G''(Ax) in L^(S^,dx). Notice that a better directionality 
of a steerable function, as measured by its auto-correlation 
function, is inevitably associated with a larger band limit A'^, 
and with a larger number T of values of m for which Gjm 
has a non-zero value for at least one value of I. Indeed, on 
the circle 5*^ as on the plane or the sphere, the smaller the 
dispersion of Ax in real space, the larger the dispersion of 
m in harmonic space. Consequently, a better directionality 
of a steerable function requires an increased number AI of 
basis functions. 

We need to understand the evolution of this direction- 
ality property of a factorized steerable function G with com- 
pact harmonic support through the kernel dilation (|2Up . The 
correlation function of two dilated versions Ga and G^' by 
factors a and a' in is defined through the scalar prod- 
uct G^^/(Ax) = {Gx,a\G^i ,a'). We consider again the case 
where the directionality coefficients are independent of I for 
I ^ N — 1 and where the azimuthal band limit for steer- 
ability is lower than the lower bound of the compact har- 
monic support of each of the two dilated versions of G: 
N ~ 1 ^ [a-^a-^Bj -Hi and AT - 1 < [a'-'^a-'^B\ + 1. 
In that limit, the correlation G£,/ (Ax) reads as 

CZ' (Ax) ^ {Ga\G,,} J2 (32) 

and appears to be simply proportional to the auto- 
correlation function of G. For a — a , this result states that 
the auto-correlations of G and Ga are proportional. The 
control of directionality of G through the auto-correlation 
function is therefore preserved through kernel dilation. For 
a 7^ a', this result essentially ensures that the kernel dilation 
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does not introduce any unexpected distortion in the shape 
of the function in real space on S^. 

In addition to the auto-correlation function, symmetry 
properties may also be imposed on the spherical harmonic 
coefficients Gim, which translate into simple directionality 
properties in real space for G. Firstly, in the framework of 
a wavelet analysis on the sphere, one generally imposes the 
symmetry relation G*„ = ( — l)™G';(_m) in order to restrict 
to real analysis functions: G{6, tp) £ R. Secondly, the con- 
straint that Gim has non-zero values only for even or odd 
values m £ Tjv, for any azimuthal band limit A'", implies that 
the function G is respectively even or odd under a rotation 
around itself by x = ^r: G{d,ip + tt) = {~l)'^-^G{e,ip). 
Thirdly, for such functions, the additional constraint that 
the spherical harmonic coefficients Gim are real for even val- 
ues of A'' — 1, and purely imaginary for odd values of A'^ — 1, 
implies that the function G is respectively even or odd un- 
der a change of sign on ip: G{0,—(p) = { — l)^~^G{6,(p). 
These symmetries are defined up to a rotation of the func- 
tion around itself by any angle x £ [0, 2tt), which amounts to 
a multiplication of the spherical harmonic coefficients Gim 
by a complex phase e~™^. The three properties discussed 
are obviously preserved through kernel dilation of factorized 
functions. They indeed only concern the directionality coef- 
ficients Si^, which are not affected by the kernel dilation. 

In summary, the auto-correlation function and addi- 
tional symmetries define directionality properties of a func- 
tion. We have shown that the directionality properties stud- 
ied are essentially preserved through kernel dilation of fac- 
torized steerable functions with compact harmonic support. 
Again, the verification of more precise directionality prop- 
erties in real space for a function designed from its spherical 
harmonic coefficients unavoidably requires a numerical eval- 
uation of sampled values of that function. 



4 WAVELETS FROM KERNEL DILATION 

In this section we begin with the derivation of a new contin- 
uous wavelet formalism from the kernel dilation with contin- 
uous scales, and for factorized steerable wavelets with com- 
pact harmonic support. We then derive the scale discretized 
wavelet formalism from the continuous wavelet formalism. 
The transition is performed through an integration of the 
dilation factor by slices. We emphasize the practical acces- 
sibility of an exact reconstruction of band-limited signals 
from a finite number of analysis scales. We also illustrate 
these developments through the explicit design of an exam- 
ple scale discretized wavelet. We finally recast the scale dis- 
cretized wavelet formalism developed in a generic invertible 
filter bank perspective. 

4.1 Continuous wavelets 

We simply consider the continuous wavelet formalism ex- 
posed in Section [5J and particularize it to the kernel dilation 
defined in Section |3] Hence the scales of analysis are still 
continuous. The translations by loq £ S'^ and proper rota- 
tions by X £ [0, 2ir) of the wavelets are still defined through 
the continuous three-dimensional rotations from relation ((2]) 
and (O. 



For application of the kernel dilation, we consider con- 
tinuous factorized steerable functions '3/ £ L^(S^,dSl) with 
compact harmonic support: 

^im = K^{l)SL, (33) 

for a continuous kernel defined by a positive real function 
K<s,{k) £ R+ and a directional split defined by the direc- 
tionality coefficients Sfm- The compact harmonic support of 
the wavelet in the interval I £ ([a^^Bj ,B) is ensured by 
a kernel K^ik) with compact support in the interval k £ 
{a~^B,B), with a compactness c{a) = a/(a — 1) £ [l,cxj): 

K^{k)^0 for k(^{a'^B,B). (34) 

The steerability of a wavelet with an azimuthal band limit 
TV in ensured by the directional split: 

•S*;^ = for all l,m with |m| ^ A'^, (35) 

with 

E = 1' (36) 

m|^min(JV-l,i) 

for all I £ N". Continuous axisymmetric wavelets A{d) with 
compact harmonic support are simply obtained by the trivial 
directional split with = SmO for all Z £ N". 

The analysis of a signal F £ L^(S^,df2) with the anal- 
ysis function ^ gives the wavelet coefficients (p, a) at 
each continuous scale a, around each point luq, and in each 
orientation x, through the directional correlation Q. The 
reconstruction of F from its wavelet coefficients results from 
relation (O. The zero mean condition (|15|) for the admissi- 
bility of * implies Kl (0) = 0. One can also set arbitrarily 
Sqo = 0. The admissibility condition (I16|l summarizes to: 

< C* = / TT^i ('^O < (37) 

J(a-^B,B) ^ 

which actually also encompasses the zero mean condition. 
The coefficients entering the reconstruction formula are 

= 8n'^Ci/{2l + 1) for I £ N". In other words, the kernel 
must formally be identified with the Fourier transform of an 
axisymmetric wavelet on the plane. 

Notice that for a factorized wavelet the directional 
correlation defining the analysis of a signal may also be un- 
derstood as a double correlation, by the kernel and the direc- 
tional split successively. The standard correlation (|9]) of the 
signal F and the axisymmetric wavelets defined by the kernel 
of vJ*, provides intermediate wavelet coefficients WJ^^{llio, a) 
on S^ at each scale a £ R!^. The spherical harmonic trans- 
form of these coefficients reads as: 

At each scale a, the directional correlation of the interme- 
diate signal obtained at that scale W^^{Loo,a) and a direc- 
tional wavelet deffned by the directional split of provides 
the final wavelet coefficients on SO (3): 

(39) 

This reasoning obviously holds independently of the steer- 
ability or compact harmonic support properties of 
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In conclusion, the definition of the kernel dilation pro- 
vides a new continuous wavelet formalism, where scales, 
translations, and proper rotations of the wavelets are all 
continuous. As the previously developed continuous wavelet 
formalism based on the stereographic dilation, it finds ap- 
plication in the identification of local directional features of 
signals on the sphere. The wavelets defined bear new proper- 
ties of compact harmonic support and steerability, which are 
preserved through kernel dilation. These properties can give 
a new insight for the analysis of local directional features. 
However, as already discussed the continuous scales required 
for the analysis prevent in practice the exact reconstruction 
of the signals analyzed from their wavelet coefficients. 

4.2 Scale discretized wavelets 

Scale discretized wavelets F can simply be obtained from 
continuous wavelets through an integration by slices of the 
dilation factor a £ R!j_. Through this transition proce- 
dure, scale discretized wavelets remain factorized steerable 
functions with compact harmonic support, and are dilated 
through kernel dilation. 

We consider the analysis of a signal F £ L^(S^,df2) 
with band limit B. The original continuous wavelet 'i' £ 
L^(S'^, dn) with a compact support is defined in the interval 
k £ {a~^B,B). The value a > 1 regulates the compactness 
c{a) of 4'. It is also taken as a basis dilation factor. The 
discrete dilation factors for the scale discretized wavelet will 
correspond to integer powers a-', for analysis depths j £ N. 

The scale discretized wavelet F £ L^(S^,df2) is thus 
defined in factorized form: 

fim^Kr{l)SL, (40) 

for a scale discretized kernel defined by a positive real func- 
tion Kr (k) £ R+ and a directional split defined by the direc- 
tionality coefficients Sf^. The directional split of F is iden- 
tified with the split of 

Sim = Sim, (41) 

also giving 

S'L =0 for all l,m with lm| ^ iV, (42) 

and 

^ jSLl' = 1, (43) 

|m|sJmin(JV-l,i) 

for Z £ N°, while Slo = 0. The exact same steerability 
properties are therefore obviously shared by the continu- 
ous wavelet and the scale discretized wavelet, independently 
of any dilation factor. The scale discretized kernel Krik) 
is obtained from the continuous kernel K\s,{k) through an 
integration by slices of the dilation factor a £ of the 
continuous wavelet formalism. 

As a first step, a positive real scaling function $r(fc) £ 
R+ of a continuous variable k £ R+, is defined which gathers 
the largest dilation factors a £ (l,oo), or correspondingly 
the lowest values of k. This generating function reads for 
A: £ R+ as: 

*r(fc) = ^ r^kliak) 

= kl ^ ^^'Hk'), (44) 



and continuously continuated at = by $r(fc) ~ 1- 
The scaling function <I>p (fc) therefore decreases continuously 
from unity down to zero in the interval k £ {a~^B, B): 

<i>r (fc) = 1 for ^ fc ^ a"^B, 

l>r (fc) £ (0, 1) for a'^ B < k < B, 

l>r (fc) = for B. (45) 

Notice that similar procedures of scale integration by slices 
were already proposed in the developmen t of corre s pond- 
ing formalisms on the plan e fPuval-Do stin et al.l Il993l: 
'Mu schietti fc Torresanil 1 19951 : IVanderghevnst fc Gobbera 
|200C 

As a second step, a si mple Littlewood-Paley decompo- 
sition (jPrazier et al.lll99lh is used to define the scale dis- 
cretized kernel Kr{k) by subtracting the scaling function 
'l>r(fc) to its contracted version 3>r(Q~^fc). This implicitly 
sets the value a as the basis dilation factor. The scale dis- 
cretized kernel also reads as an integration of the continuous 
kernel over a slice a £ (a~^,l) for the dilation factor, or 
equivalently over a slice k £ {a~^B,B) n {a~^k,k) of the 
compact support interval: 

kr (k) = $r {a~^k) - <lr (k) 

^■f J (a-^B,B)n(a-^k,k) 

The scale discretized kernel therefore has a compact support 
in the interval k £ {a~^B,aB): 

KT{k)=Q for ki{a~^B,aB). (47) 

This support is wider than for the original continuous ker- 
nel and the scaling function. The corresponding compactness 
reads as c{o?) = o?/{o? — 1) £ [l,oo). The compact har- 
monic support of the scale discretized wavelet F itself is thus 
defined in the interval I £ ([a~ B\ , [q_B]). The kernel also 
satisfies K^{Q) = 0, leading to a scale discretized wavelet F 
with a zero mean on the sphere: 

— / AQ.V{u)= 0. (48) 
47r 

The dilations by q-' of the scale discretized wavelet ob- 
tained are defined by the kernels KricPk) for any analysis 
depth j £ N. Each kernel has a compact support in the 
interval k £ {a"^^'^^'' B,cS^~^'' B) and exhibits a maximum 
in fc = B, with Kr ^ {a~-' B) — 1. The scale discretized 
wavelet T^j at each analysis depth j thus has a compact har- 
monic support in the interval I £ ( 

The property K^{0) — still ensures that each scale dis- 
cretized wavelet has a zero mean on the sphere. Notice that 
for j ^ 1, one gets a dilation factor strictly greater than 
unity a-' > 1, and the scale discretized wavelet has a band 
limit lower or equal to the assumed band limit B for the sig- 
nal F to be analyzed. At j = 0, only the values of the kernel 
in the interval / £ ([a~^-Bj , B) are of interest, as higher fre- 
quencies / are truncated by the signal F itself through the di- 
rectional correlation. One can equivalently consider that the 
compact support of the kernel is restricted to fc £ {a~^B, B) 
in the definition of the scale discretized wavelet at this first 
analysis depth j = 0. For j ^ —1, the lower bound of the 
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compact harmonic support of the scale discretized wavelet is 
larger than the band limit B. The scale discretized wavelets 
with negative analysis depths can therefore be discarded, as 
the result of their directional correlation with the signal F 
would be identically zero. 

The admissibility condition (|37p for continuous wavelets 
simply turns into a resolution of the identity below the band 
limit by a set of dilated wavelets at various analysis depths 
j, with ^ j ^ J, and a dilated scaling function at some 
total analysis depth J G N. One gets in particular for ^ 
k = l<B: 

J 

*r (a'lj + ^Ki = 1. (49) 

j=0 

The scaling function values $r(a''0 are equal to unity in 
the interval / G [0, then decrease in the in- 

terval I G ( \a "^S]), and are equal to zero 

for I ^ \a--^B'\. The kernel values KriaU) are non- 
zero only in the compact harmonic support interval I G 
( l^a-'^+^^sJ , ["a'^-^^B"! ). The scaling function typically re- 
tains the low frequency part of the signal, which will not 
be analyzed. All signal information at frequencies I ^ 

a'^^^'^'fij is kept only in the scaling function, equal to 
unity. The wavelets are equal to zero at these frequencies. 
All signal information at frequencies / ^ \a~'' B~\ is fully 
analyzed by the wavelets, while the scaling function is equal 
to zero. Intermediate frequencies are also analyzed by the 
wavelets but the scaling function is required for the recon- 
struction of the corresponding signal information. 

Let us define the maximum analysis depth Jb (a) as the 
lowest integer value such that a~'^^'^°^B ^ 1: 

Jb (a) = [log, B] . (50) 

In a case where the total analysis depth would be chosen 
strictly above Js(a), all wavelets at analysis depths j with 
J ^ j ^ Jb{c() + 1 would be identically null as their kernel 
have a compact support strictly included in the interval k G 
(0,1). The total analysis depth is consequently naturally 
limited by J ^ Js(a). In the case J = JB{a), the dilated 
scaling function evaluated at a''^^"H has a non-zero value 
only at Z = 0, 4>p(q'^^*"^Z) — 5;o, while all wavelets are 
equal to zero at Z = as they have a zero mean. Hence, the 
identity can be resolved with Jb (a) -|- 1 dilated wavelets and 
a trivial scaling function which simply retains the spherical 
harmonic coefficient _Foo out of the analysis, or equivalently 
the mean of the signal over the sphere. One gets in particular 
for < fc = Z < B: 

Sio + ^ (a'l) = 1. (51) 

4.3 Analysis and exact reconstruction 

Following the scale discretization defining the wavelets F G 
L'^(S^, df2), a new scale discretized wavelet formalism is pro- 
vided for the analysis and the exact reconstruction of band- 
limited signals. 

The analysis of a band-limited signal F G L^(S^,dr2) 
with band limit B, with a scale discretized wavelet F is 



performed by directional correlations just as in the contin- 
uous wavelet formalism. The translations by tJo G and 
proper rotations by x £ [0, Stt) of the wavelets are still 
defined through the continuous three-dimensional rotations 
from relation ([2]) and ([3|. At each analysis depth j with 
^ j ^ J ^ Jb (a) , the analysis is performed by directional 
correlations of F with the analysis functions T^j dilated 
through the kernel dilation by dilation factors a-' : 

Wr^(p,Q^) = (F,,,,|i^). (52) 

At each discrete scale a-* , the wavelet coefficients Wf{p, a-') 
still identify a square-integrable function on S0(3), and 
characterize the signal around each point wo, and in each ori- 
entation X- Once more, the direct Wigner _D-function trans- 
form of the wavelet coefficients is given as the pointwise 
product of the spherical harmonic coefficients of the signal 
and the wavelet: 

Ovf)L («0 = aTTT^'"^""- ^^^^ 

Again, the factorization relation (|40|) allows one to under- 
stand the directional correlation (|53p as a double correlation, 
by the kernel and the directional split successively. 

The reconstruction of the band-limited signal F from 
its wavelet coefficients reads in terms of a summation on a 
finite number J -f 1 of discrete dilation factors: 

F{lj) = [<i>,./F] (^) + 

J2 [ dp Wf [p, a') [r ip) L^F,,] {lo) . 

j—O SO{3) 

(54) 

The approximation F]{lj) accounts for the part of 
the signal retained in the scaling function "l>r(Q''0- ^ 
very similar way to the part of the signal analyzed by the 
wavelets, it can be written as: 

(lj) = 2tv j ^ d^o Wi {uo, [r (^o) L"^^.,]^ {uj) , 

(55) 

with W,|'(tJo, Q'') ~ {^ujo,a''W)^ and for an axisym- 
metric function >I> G L^(S^,dr2) defined by ('l?r);,„ = 
'^T{l)5mO- In the particular case where J = Jb{q), one 
gets ($r);^ = 5io5rnO and the approximation simply reduces 
to the mean of the signal over the sphere: ['^^Jb{<^)P] ~ 
(47r)~^ Jg2 dQ.F{ui). The zero mean signal is completely an- 
alyzed by the scale discretized wavelets. The operator L'^ in 
L'^(S'^, df2) in the present scale discretized wavelet formalism 
is defined by the following action on the spherical harmonic 
coefficients of functions: L'^Gim = {21 + l)G';m/87r^. This 
operator defining the scale discretized wavelets L'^T^j used 
for reconstruction is independent of F, contrarily to the op- 
erator Lif for continuous wavelets. This simply comes from 
the fact that the scale discretized wavelets are, through their 
definition (|46|). normalized by C*. 

Just as in the continuous wavelet formalism where the 
admissibility condition (|37|l is required, the present recon- 
struction formula holds if and only if the scale discretized 
wavelet satisfies the constraints (|43p , and (|49p or (|5ip . These 
constraints are automatically satisfied by construction of the 
scale discretized wavelets through the integration by slices. 
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Again, this corresponds to the requirement that the wavelet 
family as a whole, including the scaling function, preserves 
the signal information at each frequency Z £ N. 

Let us emphasize the fact that a finite number of dis- 
crete dilation factors is required for the analysis and recon- 
struction of a band-limited signal. Contrarily to the case 
of the continuous dilation factors, this allows exact recon- 
struction of band- limited signals from relation (|54|) . The 
translations and proper rotations of the wavelets are still 
defined through the continuous three-dimensional rotations. 
As discussed in Subsection 12.11 the exact reconstruction is 
achieved only for suitable pixelizations of p = (</5o,^o,x) 
which provide an exact quadrature rule for the numerical 
integration of band- limited functions on SO (3). In the case 
of non band-limited signals, an infinite number of negative 
analysis depths j ^ —1 should be added for a complete 
analysis. This would break the possibility of exact recon- 
struction. But in any case, no exact quadrature rule exists 
on S0(3) for the numerical integration of non band-limited 
functions, which already prevents an exact numerical anal- 
ysis. 

Let us also remark that scale discretized axisymmet- 
ric wavelets with compact harmonic support and dilated 
through kernel dila tion were recently introduced under the 
name of needlets |Baldi et al. I l2006l : iGuilloux et al. I l2007l : 
iMarinucci et al.ll2007^ . It is possible to show that the needlet 
coefHcients of a wide class of random signals on the sphere 
are uncorrelated in the asymptotic limit of small scales, at 
any fixed angular distance on S^. The scale discretized steer- 
able wavelets with compact harmonic support, thanks to 
their factorized form and to the choice of the kernel dila- 
tion, are also good candidates for a directional extension of 
needlets. 



4.4 Example wavelet design 

As an illustration of the transition between the continuous 
and scale discretized formalisms, we explicitly design a real 
scale discretized factorized steerable wavelet F with compact 
harmonic support, from a real continuous wavelet ^P. We 
firstly define the directional split and kernel with generic 
values for the band limit B, for the basis dilation factor a > 
1, as well as for the azimuthal band limit N of steerability. 
We then illustrate the definition for particular values. 

The directionality coefficients of 'l/ and T are identi- 
cal by the definition (|41|l . The steerability relation (|42p is 
imposed with an azimuthal band limit A''. The function is 
imposed to be real and to be even or odd both under ro- 
tation around itself by n and under a change of sign on ip. 
As discussed in Subsection l3.ll this corresponds to the con- 
straints that only the T = N values m £ Tjv are allowed, 
with Sf^ = {~^)"^ Sff^-m) 1 ™d Sfm is real for even values of 
N — 1, and purely imaginary for odd values of A^ — 1. One 
has 5oo = 0, and only the values Sf^ with 1 ^ I < B and 
^ m ^ I and m £ Tjv need to be defined explicitly. These 
values are set in order to ensure a precise structure of the 
auto-correlation function (1291) under the constraint (1431): 



N,m) 



7(iV,I) 
2 



1/2 



(56) 



values of A - 1, /3(jv,^) = [1 - (-l)^+'"]/2, and 7(iv,o = 
min(A-l, «-[H-(-l)^+']/2). The auto-correlation function 
follows as 



;G{[a-lsJ ,B) 



(57) 



When A — 1 > [a^^Sj -|- 1, the peakedness of the auto- 
correlation is generically defined by the powers 7{jv,i) of 
cos(Ax) at each value of I. This expresses the simple fact 
that the azimuthal frequency index m must always remain 
bounded in absolute value by the overall frequency index: 
m| ^ I. But the values 7(jv_;) ensure that the directionality 
coefficients are independent of I for I ^ N — 1. Hence, when 
A — 1 ^ [a~^_Bj -I- 1, the auto-correlation function takes the 
form 



C^(Ax) = ||r||^cos(^-i' (Ax), 



(58) 



with ||r|p — 5I^;g(|^(j-i3j 3) Ky{1). Its peakedness increases 
as the power A — 1 of cos(Ax). The cost for the correspond- 
ing increase of directionality with A is of course that a larger 
number of basis functions is required to steer the wavelet. 

Let us emphasize the importance of the structure (|57|) 
in the global scheme of the scale discretized wavelet for- 
malism. The azimuthal band limit A might be considered 
much smaller than the lower bound of the compact har- 
monic support interval for the first analysis depth j = 0: 
A <C Lq~^-BJ. However at each analysis depth j ^ 1, 
the compact harmonic support is defined in the interval 

I £ ([^a-^i+^'fij , ["o'l-^'fij). Hence the structure ^ of 
the auto-correlation function breaks down to ()57p at a given 
analysis depth jjv, defined as the lowest integer such that 
A — 1 > |^a~'^^-"'-'i3j -(-1. If one wants to preserve the struc- 
ture (|58[) for all dilated wavelets, the resolution of the iden- 
tity (|49|) can be used up to a total analysis depth J = j'jv — 1- 
The continuous kernel is defined from a Schwartz func- 
tion with compact support in the interval ( — 1, 1) on R as: 



K^, (k) 
(k) 



cxp 



1 - (fc) 
for t(fc)^(-l,l) 



for f(fc)£(-l,l). 



for the function 



t{k) 



^ ak~ B 
'{a-l)B 



(59) 



(60) 



with rjN = ^ for even values ofA — 1, ?7jv = i for odd 



which linearly maps the compact support interval k £ 
{a~^B,B) onto t £ (—1, 1). The function _R'*(A;) is infinitely 
differentiable for k £ R+. It notably exhibits a maximum 
at the center t{k) = of the support interval and smoothly 
drops down to zero at the interval bounds. Let us recall 
that the kernel (|59p is by definition taken as a positive func- 
tion. An overall change of sign would simply flip the sign 
of the wavelet at each point in real space. The scaling func- 
tion $r(fc) and scale discretized kernel Kr{k) follow from 
relations (|44p and (|46l) respectively. The scaling function 
for values in the interval k £ {a~^^'^-'^ B,a~-' B) for each 
analysis depth j can be obtained by numerical integration. 
Notice that the exactness of reconstruction provided by the 
formalism is not affected by such a numerical integration, 
as long as the scale discretized kernels are simply defined by 
differences of scaling functions through relation (|46p . Corre- 
sponding graphs are reported in Figure [T] for a band limit 
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Figure 1. Graphs of the continuous kernel defined in 1591 1 and 1601 1. and the corresponding scale discretized kernels obtained by differences 
of scaling functions at various analysis depths. A band limit B = 1024 and a basis dilation factor a = 2 are chosen. The continuous kernel 
Kii,{k) is represented by the continuous red line. The numerically integrated scaling function i>r{fc) is represented by the dot-dashed 
blue line. The scale discretized kernels Kii,{2^ k) are plotted as dotted black lines for the five first analysis depths j, with ^ j ^ 4. For 
j = 0, the corresponding compact support interval is cut at the band limit: k 6 (512, 1024). For 1 ^ j ^ 4 as for all larger analysis depths 
(not shown), the intervals progressively move to lower frequencies and shrink: k G (256/2'^"^', 1024/2'^"^^). At the maximum analysis 
depth j = Js(a) = 10, the compact support is shrunk to fe G (0.5, 2) and the scale discretized kernel only contains the frequency / = 1. 




Figure 2. Plots of the real scale discretized wavelets defined through relations 1561 1. 1591 1. and 1601 1. at various analysis depths. A global 
band limit B = 1024 and a basis dilation factor a = 2 are chosen, as well as an azimuthal band limit N = 3 for the steerability. Light 
and dark regions respectively correspond to positive and negative values of the functions (see value bars). The wavelets are neither 
translated, i.e. they have their central position at the North pole, nor rotated, i.e. they are in their original orientation X = (the 
meridian ip = corresponds to a vertical line passing by the North pole). The wavelets are represented at the four largest analysis 
depths, 7 ^ J ^ 10 = JB(a), identifying the four largest scales. At j = 7 (extreme left panel), j = 8 (center-left panel), and j = 9 
(center-right panel), the compact supports of the scale discretized kernels respectively contain the frequencies i = 5 to Z = 15 with a 
kernel maximum at Z = 8, Z = 3 to Z = 7 with a kernel maximum at Z = 4, and Z = 2 to Z = 3 with a kernel maximum at Z = 2. At j = 10 
(extreme right panel), the scale discretized kernel only contains the frequency Z = 1. In real space, the dispersion of angular distances 
around the central position on the sphere increases with the analysis depth, in complete coherence with the constraint (|28| l. For the 
depths j with 7 ^ j ^ 9, the lowest frequencies I are greater or equal to A'^ — 1 = 2 and the azimuthal frequency indices contained in the 
directional split are m £ {—2, 0, 2}. These wavelets all have the same directionality property as measured by an auto-correlation function 
evolving as cos^(Ax). For the depth j = 10, the scale discretized wavelet is a pure dipole (Z = 1). The azimuthal frequency index is 
restricted to m = 0, and the wavelet is simply axisymmetric with a constant auto-correlation function. 



B = 1024 and a basis dilation factor a = 2 associated with 
a standard dyadic decomposition of scales. 



Plots of the scale discretized wavelet are reported at 
various analysis depths in Figure [2] for B — 1024, a — 2, 
and for an azimuthal band limit A'' = 3 for the steerability. 
These plots notably illustrate localization and directionality 
properties of the wavelet. 



4.5 Invertible filter bank 

A scale discretized wavelet formalism with relations (|49p and 
(|51|) for factorized steerable wavelets with compact harmonic 
support can be developed by simply relying on a Littlewood- 
Paley decomposition, without any contact with the contin- 
uous wavelet formalism. One simply needs to choose any 
arbitrary scaling function satisfying relation l|45p and define 
the corresponding scale discretized kernels by differences of 
scaling functions at successive scales. 

Such invertible filter banks based on the harmonic di- 
lation were already developed in the case of axisymmetric 
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wavelets jStarck et al.|[2006lj ). and our definition of factor- 
ized steerable wavelets with compact harmonic support al- 
lows a straightforward generalization to directional wavelets 
with the kernel dilation. Also notice that the constraints 
of steerability and compact harmonic support for the scale 
discretized wavelets can technically be relaxed without af- 
fecting the Littlewood-Paley decomposition. However both 
properties are essential for the control of localization and 
directionality properties through kernel dilation. Moreover, 
in the absence of compact harmonic support, the relation 
(|49p turns into a resolution of the contracted scaling func- 
tion $p(a~^Z) which differs from unity below the band limit. 
In other words, the filter bank developed in such a case an- 
alyzes the part of the signal corresponding to its standard 
correlation with the contracted scaling function, rather than 
the signal itself. In the absence of compact harmonic sup- 
port and steerability, essential multi-resolution properties 
are also lost (see Subsection lS.ip . The memory and compu- 
tation time requirements of the algorithm for the analysis 
and reconstruction of signals therefore increase significantly 
and may rapidly become overwhelming. 

Invertible filter banks based on t he stereographic dila- 
tion have also recently been proposed (|Yeo et al.ll2006l ). but 
they do not share these essential multi-resolution properties. 



5 EXACT MULTI-RESOLUTION ALGORITHM 

In this section we identify the multi-resolution properties of 
the scale discretized wavelet formalism developed. We de- 
scribe a corresponding algorithm for the analysis and exact 
reconstruction of band-limited signals. We discuss in detail 
the memory and computation time requirements of the algo- 
rithm. Finally, an implementation of the algorithm is tested. 



5.1 Multi-resolution 

We consider the analysis and exact reconstruction of a 
band-limited signal F G L'^(S'^,dn) with a scale discretized 
wavelet F £ L^(S^, dfi), which is a factorized steerable func- 
tion with compact harmonic support. We consider a band 
limit B and a basis dilation factor a > 1. 

The signal is identified by 0{B^) spherical harmonic 
coefficients Fim- Equivalently, sampled values F{u>i) of the 
signal on a number O(B^) of points coi are generally re- 
quired in order to describe it completely. The integer i sim- 
ply indexes the points of the chosen pixelization. Notably 
exact quadrature rules for integration of band-limited sig- 
nals on S'^ with band limit B exist on equi-angular and 
Gauss-Legendre pixelizations on O(B^) points. The quadra- 
ture rules on HEALPix pixelizations on O(B^) points are 
non-exact but can be made very precise llDriscoU fc Healvl 
1 19941 : iDoroshkevich et al.|[2005al : iGorski et al.ll2005h . 

The compact harmonic support of the scale dis- 



cretized wavelet F^j is reduced in the intervals I G 
(|^a"*^+^'Bj , ["of'^'^^sj) through the kernel dilation at 
each analysis depth j. As a function on SO(3), the wavelet 
coefficients at depth j exhibit the same compact harmonic 
support as the scale discretized wavelet F^j . From relation 

(|53p, the Wigner /^-transform {Wf)^^^{a-') of the wavelet 
coefficients is indeed non-zero only in the same interval as 



the wavelet. In particular, the band limit of the wavelet 
coefficients is decreased to |"Qf'^~-'^sj at depth j. Conse- 
quently, the number of sampled values of the wavelet coef- 
ficients is reduced at each increase of the analysis depths j 
to a'^^^~-'^ X 0{B^) discrete points of the form {ujo)i(j) on 
S^, where simply indexes these points. The number of 
operations required for their computation is reduced corre- 
spondingly. Hence, the kernel dilation applied to scale dis- 
cretized wavelets with compact harmonic support provides 
a first strong multi-resolution property for the formalism. 

The steerability of the wavelet is al so important in the 
algor i thmic structure of the analysis (jWiaux et al.l |2005| . 
l2006l . I2OO7I ). beyond the fact that it ensures that direc- 
tionality properties are preserved through kernel dilation. 
Indeed, by linearity of the directional correlation (|52|l . the 
general property of steerability (|23[) is transferred from the 
wavelet to the wavelet coefficients of any signal. At each 
point (aJo)i{j) and at each analysis depth j, the wavelet co- 
efficients of a signal F with the scale discretized wavelet V^j 
are known for all continuous rotation angles x G [0, 27r) as a 
linear combination of the wavelet coefficients of F with M 
basis wavelets. As discussed, the basis wavelets can be taken 
as specific rotations T^^^^j of the wavelet on itself by rota- 
tion angles Xp ^ [0,27r), with interpolation weights given 
as simple translations by Xp of a- unique function k{x)- We 
consider wavelets for which the number of rotations required 
can by optimized to M — T ^ 2N — 1, where T is the finite 
number of values of m for which Gim has a non-zero value 
for at least one value of /. Consequently, the steerability of 
the scale discretized wavelet T^j implies a reduction of the 
number of sampled values of the wavelet coefficients to the 
T values Xp of the rotation angle, with O^p^T— 1, at 
each point (uuo) and at each analysis depth j. The num- 
ber of operations required for their computation is reduced 
correspondingly. From this perspective, steerability provides 
a second strong multi-resolution property for the formalism. 
All required sampled values of the wavelet coefficients of a 
signal with a steerable wavelet may be mapped on a sphere 
for each of the T values of Xp, at each analysis depth j. 

In summary, when multi-resolution properties of the 
formalism are fully accounted for, a reduced number of dis- 
crete points of the form p/fj) — ((ajo)i(j) , Xp) on SO(3) 
are required for the sampled values Wr^(p/(j), a-*) of the 
wavelet coefficients, where — {i{j),p} simply indexes 
these points at each analysis depth j. 



5.2 Algorithm 

The proposed algorithm works in harmonic space on and 
SO(3) in order to take advantage of the directional correla- 
tion relation (|53|) . 

Some precalculations are firstly required. The spher- 
ical harmonic coefficients (r^j)^^ of the scale discretized 
wavelets must be designed at each analysis depth j. A nu- 
merical integration can be required in order to compute the 
scaling functions ^^{a-'k) at all analysis depths from the 
spherical harmonic coefficients ^im of a continuous wavelet 
in relation (|44|l . The scale discretized kernels K^{a-' k) are 
then obtained by differences of scaling functions, and mul- 
tiplied by the directional split chosen Sf^ . 

The analysis proceeds as follows. The band-limited sig- 
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nal F is given in terms of its sampled values F{ijJi) on 
the 0{B^) discrete points coi of S^. The spherical har- 
monic coefficients Fim of the signal are computed by quadra- 
ture through a direct spherical harmonic transform. The 

direct Wigner D-function transform {W^)^^{a^) of the 
wavelet coefficients is then simply obtained by the point- 
wise product (|53p . The computation of sampled values 
{pi(j),a'') of the wavelet coefBcients requires an inverse 
Wigner D-function transform at each analysis depth j. Be- 
fore reconstruction, any suitable analysis scheme can be ap- 
plied on the wavelet coefficients, for typical purposes of de- 
noising or deconvolution. This provides altered coefficients 
Wy {pi(j),ce'')- The reconstruction proceeds trough the ex- 
act same operations as the analysis, in reverse order. The 

Wigner D-function coefficients {Wf)^^{a-') of the altered 
wavelet coefficients are computed by quadrature through a 
direct Wigner _D-function transform at each analysis depth 
j. The spherical harmonic coefficients of the reconstructed 
signal Fim are then obtained as a finite summation following 
from relations (|54p and (|55|) : 

Fim = + 

3=0 |n|!jmin(]V-l,i) 

(61) 

with 

[C^],™ = {c^'l) Fim. (62) 

In the particular case where J — Jsice), one gets trivially 
— SioSmoFoo, which corresponds to keep only the 
mean of the signal out of the analysis. 

The samples F{uJi) of the reconstructed signal are recov- 
ered by simple inverse spherical harmonic transform. If no 
alteration was applied to the wavelet coefficients, the exact 
same samples are obtained as for the original signal F. This 
exactness also obviously relies on the use of exact quadra- 
ture rules both for the direct spherical harmonic transform 
of the signal in the analysis part, and for the direct Wigner 
D-function transform of the wavelet coefficients in the re- 
construction part. This requires the choice of equi-angular 
or Gauss-Legendre pixelizations on defining the discrete 
points u!i for the sampling of the original signal, and defin- 
ing the discrete points (wo)i(j) for the sampling the wavelet 
coefficients at each analysis depth j and for each value Xp- 
Again HEALPix pixelizations provide non-exact but very 
precise quadrature rules. 

5.3 Memory requirements 

We define the storage redundancy of the algorithm as the 
ratio of the number of sampled values of the wavelet coeffi- 
cients at all analysis depths with a scale discretized wavelet, 
to the number of sampled values of the original signal itself. 
A low storage redundancy is important for achieving as low 
memory requirements as possible in a practical implemen- 
tation of the algorithm. 

The Wigner D-transform {W^)^^{a-') of the wavelet 



coefficients is non-zero only in the interval I £ 
( I^q-'^+^^bJ , ["a'^-^'sl ) at each analysis depth j. Each fre- 
quency index I is thus retained exactly twice when all analy- 
sis depths j are considered. Moreover, for a steerable wavelet 
with azimuthal band limit A'^, the index n accounting for 
the wavelet directionality in relations (|53p and (|6Hl takes 
by definition T values, with T ^ 2N — 1. On the contrary, 
the index m is only related to the signal. Consequently, the 
storage redundancy of the algorithm would be exactly 2T 
if the wavelet coefficients were to be computed in harmonic 
space only. 

However, the computation of the sampled values 
W-f {pi(^j-),a^) of the wavelet coefficients in real space on the 
discrete points p/y) = ((tJo)i(j), Xp) of SO(3) is of course 
essential for general analysis purposes. Let us recall that a 
number a'^^^~-'^ x 0{B^) of discrete points {ujo)i{j) on is 
required at each analysis depth j. This number of sampled 
values is restricted by the band limit but not by the exis- 
tence of a lower bound of the compact harmonic support. 
Thanks to the steerability, only AI — T values Xp a^re re- 
quired at each point {uio)i(j)- The storage redundancy of the 
algorithm is thus obtained by accounting for the steerability 
and summing over all analysis depths j with ^ j ^ J. In 
the most exacting case where J = Js(q), it simply reads as: 

W^B.T) (") = [l + c (l - T. (63) 

Let us recall that c(a^) — / {a^ — 1) £ [1, oo) stands for the 
compactness of the scale discretized wavelet. The number of 
the sampled values {pi{j),oi-') of the wavelet coefficients 
retained by the algorithm defines an order of magnitude of 
the memory requirements, in units corresponding to one co- 
efficient per unit of memory, as: 

[^(B.T) («) = [R^\b.t) (a) X 0{B^). (64) 

For completeness, let us emphasize that this value accounts 
for the memory requirements associated with the storage of 
the wavelet coefficients only. Memory is also required for the 
storage of the 0{B^) sampled values F{ui) of the original 
signal and the T x 0{B) values of the spherical harmonic 
coefficients 9im of the continuous wavelet, or equivalently 
of the spherical harmonic coefficients (r^j)j^ of the scale 
discretized wavelets at all analysis depths j. Additional tem- 
porary memory allocations are also necessary which depend 
on the precise implementation of the algorithm. Hence, the 
value [A'Is](B,T) (ct) is to be considered as a lower bound but 
still fixes an order of magnitude for the memory require- 
ments of the algorithm. 

In the extreme case of a large basis dilation factor 
a ^ B, the compact harmonic support of the scale dis- 
cretized wavelet essentially gets as large as the band limit, 
with a compactness c(q^) B^/{B^ — 1). The maximum 
analysis depth goes to unity, JB{ct) = 1, which implies that 
only two scales are required for the analysis. The storage 
redundancy reaches its lowest value [Rs\(b.t){o) = 2r. As 
soon as Q < B, the harmonic support of the scale discretized 
wavelet obviously gets more compact and more scales are re- 
quired. For values of the basis dilation factor very close to 
unity, the compactness gets very high and may prevent the 
wavelets to be localized enough in real space at the smallest 
analysis scale. A typically absurd value a ^ [B/{B — 1)]'^''^ 
gives c(a^) ^ B, which corresponds to a compact harmonic 



© 2007 RAS, MNRAS 000.[Tlf22l 



16 Wiaux et al. 



support selecting at maximum one frequency at a time. This 
simply reminds us of the fact that too high compactnesses 
are prohibited in the framework of a wavelet analysis. Let us 
fix ideas on practical intermediate values of a < B. Notice 
that the storage redundancy increases with the band limit 
B, as Js(a) defined in (|50l) obviously increases with B for 
a fixed value of a. We give the upper bounds in the limit 
B ^ oo and Js(a) ^ oo. A dyadic decomposition of the 
scales a — 2 corresponds to a compactness c(4) =4/3, and 
the storage redundancy is bounded by [Rs](B,T)i'2) ^ 7T/3 
for any band limit B. A steerability relation with T — 3 
hence gives a bound [Rs](^B,T)i'2) ^7. A more compact sup- 
port of the scale discretized wavelets set by a = 1.1 corre- 
sponds to a compactness c(1.21) ~ 6, and the bound on the 
redundancy rises to [7?s](b,t) (1-1) J; 7T. A value T = 3 then 
already gives [i?s](B_T) (1-1) ^ 21. 



5.4 Computation time requirements 

We define the computation redundancy of the algorithm as 
the ratio of the number of operations required for the analy- 
sis and reconstruction of a signal at all analysis depths with 
scale discretized wavelets, to the corresponding number of 
operations at the first depth (j = 0) and per azimuthal 
frequency (T = 1, as for an axisymmetric wavelet which 
contains only m = 0) . A low computation redundancy is es- 
sential for achieving as low computation time requirements 
as possible. 

The precalculation consists of the computation of the 
spherical harmonic coefficients (T^j of the scale dis- 
cretized wavelets from a continuous wavelet. At each anal- 
ysis depth j, the computation of the scale discretized ker- 
nel requires a one-dimensional numerical integration of rela- 
tion (|44l) . The corresponding number of operations required 
is independent of a as each real value k £ [1,B) is cov- 
ered exactly once by the continuous wavelets at all analysis 
depths j, whose kernels have compact supports in the inter- 
vals k G {a~^^^^^ B,a~^ B). The pointwise product between 
the scale discretized kernel and the directional split in (|40p 
requires T x 0{B) operations. As it clearly appears in the 
following, the cost of these operations is negligible relative 
to the cost of the analysis and reconstruction themselves. 
Moreover, it must be performed only once for all signals to 
be analyzed. 

The analysis at a single analysis depth j consists in a 
simple directional correlation of F with F^j on S^, leading 
to the wavelet coefficients on SO(3). The a priori number 
of operations for a naive quadrature in relation (I52|l is of 
order a®^^"-'' x 0{B^), which becomes rapidly unafl^ordable. 
Fast directional correlation algorithms based on relation 
(|53p and on the separation of the three variables of integra- 
tion on S0(3), were recently developed (jWandelt fc Gorski 
200l| : IWiaux et ahlbood : iMcEwen et "al]|2007al : IWiaux et al. 
20071 ). They allow the exact computation of the sampled val- 
ues VVr^(pj(j), a-*) of the wavelet coefficients at each analysis 
depth j through a number of operations at maximum of or- 



der a^'^~-''r X 0{B^). This number of operations is mainly 
driven by the Wigner D-function transform and naturally 
scales linearly with the number M = T of rotation angles 



Xp required by the steerability of the wavelet The recon- 
struction is symmetric to the analysis and therefore requires 
the same number of operations, in reverse order. The com- 
putation redundancy of the algorithm is thus obtained by 
accounting for the steerability and summing over all analy- 
sis depths j with ^ j ^ J.ln the most exacting case where 
J — Js(a), it simply reads as: 



-3JB(a) 



T. 



(65) 



The number of operations required by the algorithm de- 
fines an order of magnitude of the computation time require- 
ments, in units corresponding to one operation per unit of 
time, as: 



(66) 



In this expression, the impact of the compact harmonic 
support of the scale discretized wavelet is concentrated in 
c(a^) = a^/{a^ - 1) G [l,oo). 

Let us again fix ideas on practical intermediate values of 
a < B, and establish upper bounds in the limit B oo and 
Jb{oi) oo. a dyadic decomposition of the scales a = 2 
corresponds to a generalized compactness c(8) = 8/7, and 
the computation redundancy is bounded by [i?c](s,T)(2) = 
15r/7 for any band limit B. A steerability relation with T = 
3 hence gives a bound [i?c](s,T) (2) ^5 45/7 ~ 6.5. A more 
compact support of the scale discretized wavelets set by a = 
1.1 corresponds to a generalized compactness c(1.331) ~ 4, 
and the bound on the redundancy rises to [7?c](b,t)(1-1) 
5r. A value T = 3 then gives [i?c](B,T) (1-1) < 15. 



5.5 Implementation 

The proposed algorithm was implemented and tested on a 
2.2 GHz Intel Core 2 Duo CPU with 2 Gigabytes of RAM. As 
already emphasized, the choice of the pixelization on which 
the original signal F is sampled is essential to ensure the 
exactness or high precision of the mere computation of its 
spherical harmonic coefficients and hence of the whole anal- 
ysis and reconstruction process. The exactness of the pro- 
posed algorithm is simply tested by considering that the 
analysis starts at the level of the spherical harmonic coef- 
ficients Fim of the original signal, and ends at the level of 
the spherical harmonic coefficients Fim of the reconstructed 
signal F. We also tested the memory and computation time 
requirements. Let us recall that the corresponding contribu- 
tions associated with the removed direct spherical harmonic 
transform of the original signal F and inverse spherical har- 
monic transform leading to the reconstructed signal F are 
overwhelmed by the inverse and direct Wigner D-functions 
transforms required at each analysis depth j. 

Band limits B G {64, 128, 256, 512, 1024} are considered 
and the basis dilation factor is set to a = 2, hence defining a 
typical dyadic decomposition of scales. At each band limit. 



^ For the spherical harmonic transform of th e signal, fast al- 
gorit h ms exist on equi- angular pixelizations llDriscoll &: Healvl 



11994 iHealv et al.l |2003|. I2OO4) . as well as on Gauss-Legendre 
llDoroshkevich et al.ll2005al lb[) and HEALPix (Gorski et al. 200j) 
pixelizations. The corresponding number of operations required 
is at maximum of order a^'^~^' X 0{B^) thanks to the separation 
of the two variables of integration on S^. 
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B 64 128 256 512 1024 

fj. (MB) 1.3 5.3 21 84 340 

r (min) 0.019 0.092 0.73 7.0 72 

e 8.6 X lO-i'' 3.2 X lO^^^ 8.9 X lO^^^ 2.2 X lO^^^ 7.4 x 10-^2 



Table 1. Test of the implementation of the proposed algorithm for the analysis and reconstruction of signals on the sphere with scale 
discretized wavelets. Memory used /i in Megabytes (MB), as well as computation times r in minutes (min) and numerical errors e both 
averaged over five random test signals are reported, as measured on a 2.2 GHz Intel Core 2 Duo CPU with 2 Gigabytes of RAM. Five 
band limits are considered B G {64, 128, 256, 512, 1024}, and the basis dilation factor is set to a = 2. The signals are decomposed up to 
the maximum analysis depths at each band limit. The steerable wavelet used has an azimuthal band limit N = 3, with only the T = 3 
even values of the azimuthal index allowed: m G {—2, 0, 2}. 



five test signals are considered, directly defined through ran- 
dom spherical harmonic coefficients Fim with independent 
real and imaginary parts uniformly distributed in the inter- 
val ( — 1,1). The steerable wavelet F defined and illustrated 
in Subsection l4.4l is used, for an azimuthal band limit N = 3. 
It only contains the T = 3 even values m € { — 2,0,2}. As 
discussed in Subsectionl3.1l the basis functions for the steer- 



ability can be chosen as the three rotated versions F^^ with 
Xp = 7rp/3 for ^ p ^ 2, and the function k{x) follows 
accordingly. The analysis is performed up to the maximum 
analysis depth for each band limit: J64(2) = 6, Ji28(2) = 7, 
J256(2) = 8, J5i2(2) = 9, and Jio24(2) = 10. The numeri- 
cal error associated with the algorithm is evaluated as the 
maximum absolute value, for all values of / and m, of the 
difference between the original and reconstructed spherical 
harmonic coefficients: e = max;,™ \Fim — Fim\- The algo- 
rithm is coded with double precision numbers, which sets 
the unit of memory for the storage of coefficients to 8 bytes. 

The memory used fi, as well as computation times r 
and numerical errors e both averaged over the five random 
test signals are reported in Table [T] The values reported 
respectively illustrate the memory requirements (|64|) . the 
computation time requirements ((66}, and the exactness of 
reconstruction of the proposed algorithm. 



6 ASTROPHYSICAL APPLICATION 

In this section we emphasize an important astrophysical ap- 
plication of the wavelet formalism defined and implemented, 
discussing in some detail the issue of the detection of cos- 
mic strings through the denoising of full-sky CMB data. We 
firstly introduce the question of the existence of topological 
defects in the Universe. We highlight the non-Gaussianity of 
the component of the CMB signal induced by cosmic strings 
and justify a wavelet decomposition of the data as a way to 
enhance the sparsity of the wavelet coefficients of the string 
signal. We then propose a denoising method based on a sta- 
tistical model of the wavelet coefficients of the string signal. 
We also emphasize the need for a precise test allowing one 
to set a confidence level on the string signal reconstructed 
from the denoised wavelet coefficients. 



6.1 Topological defects 

Observations of the CMB and of the Large Scale Struc- 
ture (LSS) of the Universe have led to the definition of a 
concordance cosmological model. The full-sky data of the 
WMAP experiment have played a domi nant role in develop- 
ing this precise pic t ure of the Universe feennett et al.| 2003 ; 
'Spergel ot al.'2003; Hinshaw ot al.'2007; Sporgcl ct al. 2003; 
,Hinshaw ct al. 2008; Komatsu ct al.. ,2008 ). In this frame- 
work, the cosmic structures originate largely from Gaus- 
sian adiabatic perturbations seeded in the early phase of 
infiation of the Universe. However, cosmological scenarios 
motivated in the context of theories unifying the funda- 
mental interactions suggest the existence of topological de- 
fects resulting from phase transitions at the end of infia- 
tion. These defects would have participated to the formation 
of the cosmic structures. While textures are more or less 
axisy mmetric, cosm ic strings a re a line-like version of de- 
fects (IVilenkin fc Shellard, 19SM| : iHindmarsh fc KibbldHooBI : 



iTurok fc Spergellliggd ). Even though observations largely fit 
with an origin of the cosmic structures in terms adiabatic 
perturbations, room is still available for the existence of a 
small fraction of topological defects. Moreover, fundamental 
string theory predicts th e existence of cosmic s trings in the 
"brane-world" scenario IIDavis fc Kibblell2005f ). As a con- 
sequence, the issue of the existence of topological defects 
represents today a central question in cosmology. 

Textures would induce hot and cold spots in the CMB 
with typical a ngular sizes of s everal degrees on the ce- 
lestial sphere ([Turok fc Spergell [l990). A recent analysis 
(|Cruz et al.]|2007l ') of the WMAP data showed that the cold 
spot detected at {0,^p) = (147°, 209°) in Galactic spheri- 
cal coordinates, is satisfactorily described by a texture with 
an angular size around 10°. The main signature of cos- 
mic st rings in the CMB is kn own as the Kaiser-Stebbins 
effect (|Kaiser fc Stebbin3ll984h . characterized by tempera- 
ture steps along the strings, with a typical angular size be- 
low 1° on the celestial sphere. Constraints have been set 
on a possible string contribution in terms of upper lim- 
its on the so-called string tension Gn, where G stands for 
the gravitational constant. The string tension sets the over- 
all amplitude of the string contribution. These constraints 
mainly come from the analysis of the string contribution to 
the overall CM B angular power spectrum (|Contaldil Il999l : 
IWvman et aLir2005, 20od iBevis et al.ll2007h . Very few algo- 
rithms have been designed for the explicit identification of 
cosmic strings through the Kaiser-Stebbins effect on full-sky 
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data (|jeong fc Smootll2005l : ILo fc Wrightll2005l ). No strong 
detection of cosmic strings has ever been reported. 

Current CMB experiments, among which WMAP, 
achieve an angular resolution on the celestial sphere of the 
order of 10 arcminutes, corresponding to a limit frequency 
B ~ 2 X 10''. These experiments constrain a possible string 
signal to be largely dominated by the standard Gaussian 
CMB contribution at the frequencies I probed, but it might 
nevertheless become a dominant contribution at higher fre- 
quencies, due to the slow decay of the corresponding an gular 
power spectrum (,Fraisse et al.,.2007l : iBevis et al.ll2007^ . The 
Planck experiment will provide full-sky CM B data at a reso - 
lution of 5 arcminutes, i. e. with _B ~ 4 x 10 l|Bouchetll2004 ). 
Important new information relative to a cosmic string signal 
will therefore be available. 

In this perspe ctive, we sketch in t he following a new sta- 
tistical approach (jWiaux et al.|[200^ for the identification 
and reconstruction of cosmic strings through the denoising 
of full-sky CMB data. It is specifically considered in the 
framework of the scale discretized steerable wavelet formal- 
ism on the spheral. Further refinement of this approach, as 
well as its precise implementation, its application to CMB 
data, and its comparison with other detection algorithms, 
are the subjects of a future work. 

6.2 Non-Gaussian string signal 

The standard component of the CMB signal induced by 
adiabatic perturbations is a Gaussian signal on the sphere 
with a known angular power spectrum in a given cosmo- 
logical model. Topological defects, and in particular cosmic 
strings, induce a non-Gaussian component of the CMB sig- 
nal with characteristic features defined at specific positions 
and scales. The corresponding angular power spectrum ex- 
hibits a fixed characteristic shape, with a slow decay at high 
frequencies. The complete non-Gaussian statistical distribu- 
tion of a string signal and the corresponding angular power 
spectrum can indeed be de duced from simulations in the 
chose n cosmological model (^Fr aisse et al. I l2007l : Ib evis et al.l 
I2OO7I ). up to an overall amplitude of the string contribution 
as set by the unknown string tension Gfi. These two sta- 
tistically independent components simply add linearly. We 
consider the perturbations of the signals around their sta- 
tistical mean. Instrumental Gaussian white noise with zero 
mean also unavoidably adds as an independent component, 
setting the limited sensitivity of the experiment considered. 

* Notice tha t experiments such as the Arcmin ute Microkelvin Im- 
ager (AMI) ll Jones et al.ll200lj Barker et al.|[2006). the Atacama 
Cosmology Teles cope (ACT) ( Kosowsky 2006), or the South Pole 
Telescope (SPT) l lRuhl et al.12004,') will map the CMB at a resolu- 
tion around 1 arcminute, i. e. with B ~ 2 X 10^ . The corresponding 
prospects for the detection of strings are thus improved relative to 
Planck data, but these experiments will provide observations of 
small portions of the celestial sphere only. Specific algorithms for 
the identification of cosmic strings in CMB data on planar patches 
must be considered. As a formalis m of scale discretized s teerable 
wavelets also exists on the plane l lSimoncelli et al.lllQQ^) . it can 
also be used for the identification and reconstruction of cosmic 
strings through the denoising of CMB data on planar patches, 
in a statistical app roach analogous to the one proposed below 
llWiaux et al.ll2008r) . 



We leave apart any issue of deconvolution of the experi- 
mental beam and also discard problems of contamination of 
CMB data by foreground emissions. In the perspective of 
the detection of cosmic strings, the non-Gaussian compo- 
nent from strings represents the signal to be identified and 
reconstructed, while the Gaussian components can be seen 
as a statistically independent Gaussian noise. The overall 
signal F reads as the sum of the string signal and the noise 
in terms of a linear combination 

F (lo,) = asFs (lo,) + F„ (lo,) , (67) 

where Fs represents the string signal for a string tension 
Us = G/i normalized to unity, and Fn represents the noise. 
The zero mean signals F, Fs, and F„ are considered to have 
a band limit B, related to the resolution of the experiment 
under consideration, and a number 0{B^) of points uji are 
required for their precise description. 

In first approximation we can fix the cosmological pa- 
rameters at their values in the concordance cosmological 
model. This fixes the angular power spectra of the noise 
Fn and of the normalized string signal Fs- 



6.3 Sparse wavelet coefficients 

Wavelets are by construction filters with zero mean (see 
(|48|) ). As such, they generically enhance discontinuities, and 
reduce smooth patterns in the signal analyzed. The non- 
Gaussian string signal characterized by temperature steps 
typically has a sparse expansion in terms of wavelets. It in- 
deed only exhibits a small number of wavelet coefficients of 
large absolute value at the specific positions of the strings 
and at their characteristic scales. On the contrary, the Gaus- 
sian contributions are characterized by smooth patterns de- 
signed by their angular correlation functions. Their expan- 
sion in terms of wavelets is not sparse. The sparsity of the 
wavelet coefficients of the string signal justifies the wavelet 
decomposition of the data. Moreover, directional wavelets 
(i.e. with an azimuthal band limit A*' > 1) particularly ap- 
ply for an efficient detection of the localized directional fea- 
tures associated with the Kaiser-Stebbins effect. Indeed, the 
more similar the filter to the signal signatures, the better it 
magnifies these signatures. Correspondingly, the detection of 
textures would more naturally follow from an analysis with 
axisymmetric wavelets {N = 1). Finally, as emphasized al- 
ready, the scale discretization of the wavelets is essential for 
the reconstruction of the signal. In conclusion, the denoising 
procedure will be more efficient when applied to the wavelet 
coefficients of the signal observed F decomposed with a scale 
discretized steerable wavelet. 

A suitable scale discretized steerable wavelet P is thus 
chosen with a given directionality set in terms of an az- 
imuthal band limit A''. The basis dilation factor a and the 
total analysis depth J are also chosen, in order to opti- 
mize the number of analysis depths in the range of frequen- 
cies I concerned by the string signal. Let us recall that the 
sampling of the wavelet coefficients of a signal is defined 
on the points pifj-^ = {{ujo)i(j),Xp) on SO(3). The value 
^0) = simply indexes these points at each anal- 

ysis depth j with ^ j ^ J ^ Jb{oi). They correspond to 
the points (ajo)i{j) on for each value of the rotation an- 
gle Xpi with O^p^^T — 1, as required by the steerability 
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(see Subsection [STTJ . By linearity of the wavelet decomposi- 
tion (|52|) . the wavelet coefficients W-p of the overall signal F 
read, in terms of the wavelet coefficients of the normalized 
string signal Fs and of the noise Fn , as 



(68) 



with OsVKp 



The wavelet coefficients W-p , Wj 



and Wp " have zero statistical means just as the correspond- 
ing signals. 



6.4 Statistical model 

In a training phase, the statistical distributions of the 
wavelet coefficients of a pure noise Fn and of a pure nor- 
malized string signal Fs must be identified. 

For the noise F„, the wavelet coefficients remain Gaus- 
sian by linearity. Assuming the statistical isotropy of the 
noise, the probability density functions of the zero mean 
wavelet coefficients VKjf" depend on the analysis depth j 



but are independent of Pi(j)- 



exp 



(69) 



The variances (o'f")^ can be inferred from the known angu- 
lar power spectrum of the noise in the range of frequencies 
probed by the wavelets at the different analysis depths. 

For the normalized string signal Fs, a Monte Carlo anal- 
ysis based on stri ng signal simulations (|Bevis et al.l l2007l : 
iFraisse et al.|[2007l ) is required to fit a non-Gaussian model 
of the probability density function at each depth j. The 
more reliable the simulations, the better the model for the 
probability density functions. As a first approximation, the 
computation of the variance and kurtosis of the wavelet co- 
efficients allows one to fit a generalized Gaussian distribu- 
tion at each depth. Assuming the statistical isotropy of the 
string signal, the probability density functions of the zero 
mean wavelet coefficients Wp" again depend on the analysis 
depth j but are independent of Pi(j) : 



■ exp 



(70) 



The parameters Uj relate to the standard deviations " 
of the distributions. The corresponding variances (o"/")^ re- 
flect the angular power spectrum of the string signal in the 
range of frequencies probed by the wavelets at the difi'erent 
analysis depths. The parameters hj obviously measure the 
peakedness of the distributions, and relate to their kurtoses. 
At the analysis depths corresponding to the characteristic 
range of frequencies concerned by the string signal, the spar- 
sity of the wavelet coefficients can be associated with peaked 
distributions f^" with heavy tails (i.e. with kurtoses larger 
than 3, or values < hj < 2) relative to a Gaussian distribu- 
tion {i.e. with a kurtosis equal to 3, or a value ^Gaussian = 2). 



6.5 Denoising 

The identification and reconstruction of a string signal from 
real data can then be implemented as follows. 



Firstly, the overall signal F is decomposed with the cho- 
sen scale discretized steerable wavelet F, which gives the 
wavelet coefficients Wp through relation (|52[) . 

Secondly, the string tension associated with a still hy- 
pothetical string signal is estimated in relation (|67|l . A pre- 
cise approach based on the analysis of the angular power 
spectrum of the real data could be adopted, consisting in 
a likelihood analysis involving all cosmological parameters 
including the string tension. This standard approach was 
used to obtain the cu rrent constraints on the string ten- 
sion (|Bevis et al.ll2007l ). together with a reassessment of the 
other cosmological parameters. The distributions f^" and 
f^" should then also be reassessed according to the modified 
angular power spectra for F„ and Fs, respectively. Ifowever, 
in the approximation considered, all cosmological parame- 
ters are fixed at their values in the concordance cosmological 
model throughout the analysis. The angular power spectra 
of the noise Fn and of the normalized string signal Fs are 
thus kept invariant. A rough estimation of the string tension 
as — Gfi is obtained from a least squares method based on 
the variances of the wavelet coefficients of F, Fs, and Fn at 
all depths j. It is primarily intended to serve to the denoising 
itself and not as a final estimation of the string tension. The 
wavelet decomposition here simply helps to bin the values of 
the power spectra before defining the constraints. By statis- 
tical independence, the variances (cyf)^ of the overall signal 
F at each depth j read, in terms of the variances {o'f")^ of 
the normalized string signal Fs and of the variances (fj^")^ 
of the noise Fn, as 



(71) 



with a^{a^ 



asFss2^ The wavelet coefficients Wp are 



thus assumed to satisfy relation 



for an estimation Us of 



the exact value as. The statistical distributions fj " for the 
wavelet coefficients Wjf" of the noise Fn are identified as 
in (|69p . The statistical distributions f^"^^ for the wavelet 
coefficients Wp"^" of the string signal dsFs are as in (|7U|I 
with Uj — > dsUj and hj left invariant. By Bayes' theorem, 
the posterior probability distribution function jj^^^"!^' at 
each depth j for the wavelet coefficients Wp"^" given the 
observed values Wf reads as 



J ] 



(72) 



Notice that one could also easily account for a flexibility in 
the overall amplitude of the standard Gaussian component 
of the CMB, associated with the cosmological parameter erg, 
in the same least squares approach. 

Thirdly, from the identifled posterior probability, the 
wavelet coefficients Wp"^' of the string signal are estimated 
to values Wp'^' , separately at each point Pi{_,) for each anal- 
ysis depth j. For example, in a maximum a posteriori ap- 
proach, this estimation is defined as the value which max- 
imizes the posterior probability, while in a Bayesian least 
square approach, it is defined as the expectation value of 
the posterior probability. 

Finally, the estimated string signal asFs is recon- 
structed from the denoised wavelet coefficients Wr'^' 



through relations (|54|) and (|55|) . The string network im 
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printed in the analyzed CMB data is readily mapped as the 
magnitude of gradient of the reconstructed signal. 



6.6 Confidence level of detection 

For the reasons discussed above, scale discretized steerable 
wavelets should represent a very powerful tool for the identi- 
fication and the reconstruction of the string signal buried in 
the standard Gaussian component of the CMB and in instru- 
mental noise. The statistical approach considered for the de- 
noising procedure at each point pi(j) for each analysis depth 
j specifically accounts for the shape of the power spectra of 
the string signal and of the Gaussian components. It also 
accounts for the peakedness of the non-Gaussian string sig- 
nal, characterized by a large kurtosis at each analysis depth 
j characteristic of the string signal. 

After reconstruction, a hypothesis test can be set up in 
order to assess if the estimated string network indeed arises 
from a string signal with high probability. For example, the 
kurtosis of the magnitude of gradient of the reconstructed 
signal can be compared to the corresponding reconstructed 
kurtosis of combinations (|67p of a string signal with noise, 
through Monte Carlo analyses for various string tensions. 
More statistics may also be combined, such as the variances 
and kurtoses at different resolutions of the magnitude of gra- 
dient of the reconstructed signal, in order to provide a more 
robust hypothesis test. This procedure is also intended to 
provide an estimation of the string tension Ua = G/i from 
the denoised signal, more precise than the original estima- 
tion 0,3 which first served to the denoising. 

Let us finally notice that the Canny algorithm (|Cannvl 

1 19861 ) was recently proposed for the d etection o f co smi c 
strings in CMB data on planar patches (l Amsel et al. 2 0071 ). 
It consists of an edge detection in the map of the magni- 
tude of gradient of the original signal, independently of any 
denoising approach. This method could be implemented for 
the analysis of full-sky CMB data, and compared to our al- 
gorithm in order to assess their relative performances. In 
the context of our denoising approach, this edge detection 
might actually be applied to the magnitude of gradient of 
the denoised signal, in order to count the string segments 
obtained, instead of computing the corresponding kurtosis. 
The confidence level of the detection could then be assessed 
by comparison with the corresponding counts for pure noise 
through a Monte Carlo analysis. 



7 CONCLUSION 

We have derived a scale discretized wavelet formalism for the 
analysis and exact reconstruction of band-limited signals on 
the sphere with directional wavelets. The combination of the 
two properties of exact reconstruction and directionality was 
lacking in the existi ng wavelet formalisms. As f or the for- 
malism developed b y lAntoine fc Vanderghevnsil (|l999l ') and 
IWiaux et al.l l|2005h , the translations of the wavelets at any 
point on the sphere and their proper rotations are still de- 
fined through the continuous three-dimensional rotations. 
But the wavelets are factorized steerable functions with com- 
pact harmonic support, and they are dilated through a ker- 
nel dilation directly defined in harmonic space. 



As an intermediate step, a continuous wavelet formal- 
ism was obtained. This by-product of our developments can 
be understood as an alternative approach for the analysis of 
signals, with wavelets bearing new compact harmonic sup- 
port and directionality properties. However, the continuous 
range of scales required for the analysis still prevents in prac- 
tice the exact reconstruction of the signals analyzed from 
their wavelet coefficients. 

The scale discretized wavelet formalism results from an 
integration by slices of the dilation factor of the continu- 
ous formalism. It allows in practice the exact reconstruction 
of band-limited signals from their wavelet coefficients with 
a finite number of scales. It can also be derived indepen- 
dently of the continuous wavelets, and can be understood as 
a generalization of existing invertible filter bank methods. 
The multi-resolution properties of the formalism were iden- 
tified and a corresponding exact algorithm was described. 
The memory and computation time requirements were dis- 
cussed and an implementation was tested. 

This formalism is of interest in a large variety of fields, 
notably for the denoising or the deconvolution of signals on 
the sphere with a sparse expansion in wavelets. It typically 
concerns signals identified by directional features at specific 
positions and scales. In astrophysics, it finds a particular 
application for the identification of localized directional fea- 
tures in CMB data, such as the imprint of topological de- 
fects, in particular cosmic strings, and for their reconstruc- 
tion after separation from the other signal components. In 
this context, we have discussed a new statistical approach 
for the detection of cosmic strings through the denoising 
of full-sky CMB data. This application is the subject of a 
future work. 
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